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1  Executive  Summary 

We  carried  out  a  synoptic  study  of  P  coda  waves  recorded  at  the  medium- 
aperture  array  located  in  Chiang  Mai,  Thailand  (CMAR).  In  particular,  we 
focused  on  seismograms  recorded  at  distances  of  13°  —  30°  in  which  multiple 
arrivals  are  created  by  velocity  discontinuities  and  gradients  in  the  upper 
mantle.  These  arrivals  can  complicate  the  goals  of  verification  seismology 
by,  for  instance,  causing  confusion  in  the  identification  of  depth  phases.  Ar¬ 
ray  processing  techniques  allow  incoming  energy  to  be  sorted  according  to 
slowness,  and  so  arrays  are  powerful  tools  in  identifying  the  origin  of  later- 
arriving  coda  waves.  In  general,  energy  with  vector  slowness  similar  to  the 
first  arrival  primarily  comes  from  near-source  scattering  (including  depth 
phases),  energy  with  significantly  different  vector  slowness  primarily  comes 
from  near-receiver  scattering,  and  energy  arriving  along  the  great  circle  path, 
but  with  slightly  different  ray  parameter,  is  possibly  related  to  upper  mantle 
triplications. 

To  help  understand  the  capabilities  of  CMAR  in  detecting  regional  dis¬ 
tance  P  coda  waves  and  resolving  their  slowness,  we  first  made  a  thorough 
study  of  the  ambient  noise  held.  As  expected  from  the  CMAR  geometry, 
we  found  that  a  frequency  band  centered  near  1  Hz  provided  the  optimal 
separation  between  signal  and  noise  coherence.  Although  there  were  large 
variations  from  event  to  event,  in  this  band  correlation  coefficients  averaged 
0.5-0. 6  for  signals  and  0.0  for  noise.  This  corresponds  to  SNR  gains  of  about 
3.2  with  standard  linear  array  processing.  At  frequencies  below  about  1.4  Hz, 
the  ambient  noise  was  strongly  anisotropic  with  respect  to  both  apparent  ve¬ 
locity  and  backazimuth.  Noise  preferentially  arrived  with  apparent  velocities 
centered  near  4  krn/s  (higher  mode  Rayleigh  waves)  and  25  km/s  (teleseismic 
P  waves),  while  the  ring  of  slowness  space  in  between  (regional  distance  P 
waves)  was  relatively  quiet.  The  Rayleigh  noise  arrived  preferentially  from 
the  southwest  and  east,  while  the  northerly  directions  corresponding  to  the 
main  Asian  landmass  were  especially  quiet. 

Our  primary  data  set  consisted  of  over  950  seismic  events  that  occurred 
from  1994-2004  at  distances  of  13°  —  30°  from  CMAR  with  continental  paths 
in  South  Asia.  Accurate  hypocenters  for  these  events  were  obtained  from 
the  EHB  catalog  and  all  the  data  were  visually  inspected  and  confirmed  to 
be  of  high  quality.  We  devised  an  algorithm  that  performed  a  sliding  win- 
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dow  slowness  analysis,  essentially  converting  the  18  time  series  of  ground 
velocity  recorded  by  the  elements  of  CMAR  into  4  times  series  representing 
beam  power,  beam  coherence,  apparent  velocity,  and  backazimuth.  We  used 
a  second  algorithm  to  sort  through  the  four  derived  time  series  and  identify 
later  arrivals  that  met  specific  criteria.  In  this  way,  the  later  arrivals  were  se¬ 
lected  automatically  and  objectively.  We  also  binned  and  stacked  envelopes 
of  beam  power  in  order  to  determine  average  coda  shapes  as  a  function  of 
distance. 

We  found  that  nearly  all  significant  waves  arriving  within  the  first  thirty 
seconds  of  the  P  wave  traveled  along  the  great  circle  path.  This  agrees  with 
previous  array-based  studies  of  teleseismic  P  coda  that  have  generally  found 
near-source  scattering  to  be  dominant  over  near-receiver  scattering.  Near¬ 
source  scattering  was  especially  strong  at  distances  of  27°  —  30°,  for  which  it 
was  common  to  observe  long  trains  of  scattering  energy  leading  to  average 
envelopes  of  beam  power  that  were  nearly  flat  for  30  s  after  the  first  arrival. 
A  potential  explanation  for  these  observations  is  near-source  scattering  of  Rg 
to  P  energy,  perhaps  by  topography. 

We  examined  various  schemes  for  stacking  array  beams  and  in  general 
found  it  difficult  to  make  coherent  images  of  triplicated  arrivals.  We  at¬ 
tribute  this  to  (1)  the  large  amount  of  near-source  scattered  energy  that 
exists  at  the  frequencies  (~1  Hz)  for  which  CMAR  is  effective  as  an  ar¬ 
ray,  (2)  the  relatively  small  aperture  (~10  km)  and  spacing  of  CMAR  that 
limits  its  slowness  resolution,  and  (3)  regional  variations  in  transition  zone 
structure  that  cause  beam  stacks  to  become  defocused.  The  only  conclusive 
triplicated  arrivals  we  observed  were  associated  with  events  that  occurred 
14°  —  16°  away  from  CMAR  in  Tibet.  Compared  to  the  initial  P  waves 
these  arrivals  were  6-8  s  later,  3-5  times  larger,  and  had  ray  parameters  2- 
3  s/deg  smaller.  These  properties  are  consistent  with  waves  that  have  either 
dived  beneath  the  410  km-discontinuity  or  been  post-critically  reflected  from 
it.  The  differential  times  are  slightly  larger  than  predictions  from  IASP91, 
which  may  imply  a  slightly  deeper  410-km  discontinuity  beneath  eastern  In¬ 
dia  and  Burma;  however,  the  first  arrivals  have  negative  absolute  residuals 
themselves  and  so  it  may  be  that  the  upper  portion  of  the  upper  mantle  in 
this  area  is  seismically  fast. 
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2  Introduction 

2.1  Motivation 

The  current  emphasis  on  detecting  and  locating  low-yield  nuclear  explosions 
makes  understanding  regional  distance  (13°  —  30°)  seismic  data  especially  im¬ 
portant.  At  telcseismic  distances  signals  are  often  too  weak  to  be  recorded 
and  at  local  distances  there  are  often  few  (or  even  no)  accessible  seismic 
stations.  Bnt  body  wave  propagation  at  regional  distances  is  complicated 
because  the  waves  turn  in  the  upper  mantle  and  strongly  interact  with  the 
seismic  discontinuities  and  sharp  radial  gradients  in  the  mantle  transition 
zone  (MTZ).  These  interactions  result  in  the  appearance  of  several  discrete 
geometric  arrivals,  often  referred  to  as  a  triplication  or  triplicated  arrivals. 
Here  we  use  the  more  general  term  of  multipathing  to  describe  this  phe¬ 
nomenon.  Although  multipathing  is  sometimes  used  specifically  in  the  con¬ 
text  of  lateral  variations  in  velocity,  we  emphasize  that  in  this  study  all  of 
the  interpretation  and  modeling  assumes  ID  wave  propagation. 

Multipathing  complicates  the  goals  of  verification  seismology  in  several 
ways.  At  certain  distances  the  first  arrival  may  be  5-10  times  smaller  than  a 
secondary  arrival.  Depending  on  which  peak  is  chosen,  the  station  estimate 
of  nib  can  be  significantly  biased.  It  is  also  possible  that  a  small  first  arrival  is 
obscured  by  noise,  leading  an  analyst  to  pick  the  larger  secondary  arrival  as 
the  nominal  first  arrival.  The  corresponding  travel  time  anomaly  then  biases 
the  location  estimate.  The  existence  of  multiple  arrivals  can  also  lead  to  con¬ 
fusion  in  identifying  depth  phases,  which  are  one  of  the  most  effective  tools 
in  discriminating  earthquakes  from  explosions.  Furthermore,  multipathing 
could  lead  to  the  inference  of  a  complicated  source  time  function,  or  multiple 
discrete  sources,  when  in  fact  only  a  single,  simple  source  exists.  On  the 
other  hand,  if  multipathing  is  properly  recognized  the  differential  times  can 
lead  to  increased  location  accuracy  (e.g.  Ringdal,  1981). 

It  is  well-known  that  the  structure  of  the  MTZ  varies  from  region  to  re¬ 
gion.  A  recent  comprehensive  study  of  receiver  functions  found  that  MTZ 
thickness  varies  globally  by  ±  35  km  (Lawrence  &  Shearer,  2006).  Thickness 
is  defined  as  the  difference  between  the  depths  to  the  410-km  and  660-krn 
discontinuities,  and  is  more  robust  than  estimates  of  topography  on  either 
discontinuity.  In  addition  to  these  two  standard  global  discontinuities,  there 
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is  evidence  for  many  other  discontinuities  that  exist  at  least  regionally  and 
perhaps  globally.  These  include  the  Lehmann  discontinuity  near  a  depth  of 
220  km  (e.g.,  Dziewonski  &  Anderson,  1981;  Karato,  1992;  Gaherty  &  Jor¬ 
dan,  1995;  Deuss  &  Woodhouse,  2004),  the  X-discontinuity  at  depths  near 
300  km  (e.g.,  Revenaugh  &  Jordan,  1991;  Ganguly  &  Frost,  2006),  and  a 
discontinuity  near  520  km  (e.g.,  Shearer,  1990;  Ryberg  et  al.,  1997;  Deuss 
&  Woodhouse,  2001).  Furthermore,  there  is  evidence  for  extreme  regional 
variation  in  the  velocity  structure  just  above  the  410-km  discontinuity,  per¬ 
haps  indicative  of  partial  melt  (Revenaugh  &  Sipkin,  1994;  Song  et  ah,  2004). 
Therefore,  each  region  of  interest  in  verification  seismology  must  be  indepen¬ 
dently  studied  and  in  a  sense  calibrated  if  upper  mantle  multipathing  is  to 
be  properly  understood. 

In  this  project  we  target  the  continental  region  of  South  Asia  as  seen  by 
the  small  aperture  seismic  array  in  Chiang  Mai,  Thailand  (CMAR).  We  use 
an  array  station  because  it  allows  for  a  large  increase  in  signal-to-noise-ratio 
(SNR),  especially  at  the  high  frequencies  common  in  regional  seismograms 
of  small  events.  Since  CMAR  has  been  active  for  over  ten  years,  this  leads  to 
a  relatively  large  database  of  earthquake  waveforms  allowing  for  a  synoptic 
study.  An  array  station  also  provides  the  capability  of  using  slowness  analysis 
to  separate  multiply  arriving  waves  by  ray  parameter,  or  in  other  words,  to 
apply  wavenumber  Liters  to  the  data.  While  depth  phases  and  near-source 
scattered  waves  are  expected  to  have  ray  parameters  close  to  the  first  arrival, 
multipathing  created  by  MTZ  structure  is  expected  to  produce  arrivals  that 
vary  by  as  much  as  3-4  s/deg  relative  to  the  first  arrival.  We  also  note  that 
CMAR  is  located  at  regional  or  far-regional  distances  from  several  known 
test  sites,  including  Lop  Nor,  China  (24.7°);  Pokhran,  India  (26.4°);  Ras 
Koh,  Pakistan  (32.9°);  and  Chik-tong,  North  Korea  (34.4°). 

2.2  Wave  Propagation  in  the  Mantle  Transition  Zone 

Example  paths  for  rays  traversing  the  mantle  transition  zone  are  shown  in 
the  top  panel  of  Figure  7.1.  Here  we  use  the  IASP91  velocity  model  (Kennett 
&  Engdahl,  1991)  and  consider  a  range  of  ray  parameters  from  8.0  s/deg  to 
15.0  s/deg,  equally  spaced  in  increments  of  0.05  s/deg.  The  head  wave  Pn 
is  not  naturally  simulated  with  geometric  ray  tracing,  however  a  rich  variety 
of  refractions  and  reflections  are  visible.  The  focusing  of  energy  near  1°  —  2° 
is  created  by  postcritical  reflections  from  the  Moho  (PmP)  and  is  commonly 
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observed  in  seismograms  recorded  at  local  distances.  More  relevant  to  this 
study  is  the  focusing  of  energy  at  15°  —  20°  created  by  returns  from  the 
410-km  and  660-km  discontinuities.  This  is  also  the  most  extreme  region 
of  multipathing,  with  IASP91  predicting  7  distinct  arrivals  (not  including 
non-geometrical  contributions  to  Pn)  at  some  of  these  distances.  Generally, 
by  30°  the  multipathing  ends,  however  for  velocity  models  with  a  strong  dis¬ 
continuity  near  220  km,  such  as  PREM,  asthenospheric  phases  are  predicted 
out  to  nearly  40°. 

In  the  middle  panel  of  Figure  7.1  we  present  ray  paths  for  all  the  geomet¬ 
ric  P  waves  at  17°.  IASP91  predicts  5  distinct  arrivals,  two  associated  with 
the  low  velocity  asthenosphere,  one  turning  above  the  410,  one  reflecting  off 
of  the  410,  and  one  diving  beneath  the  410.  The  three  shallowest  rays  arrive 
within  1.3  s  of  one  another  and  differ  by  at  most  1  s/deg  in  ray  parameter. 
Therefore,  it  would  be  difficult  to  discriminate  among  these  arrivals,  even 
with  data  from  a  high  quality  array  station.  This  is  especially  true  because 
near-source  scattered  energy,  which  often  contributes  to  the  coda,  is  expected 
to  have  a  similar  ray  parameter.  However,  the  two  deeper  rays  arrive  at  ray 
parameters  2-3  s/deg  smaller  than  the  three  shallow  rays,  and  so  even  though 
the  predicted  times  are  somewhat  close,  a  wavenumber  filter  may  separate 
the  later  arrivals  from  coda  and/or  depth  phases  of  the  first  group  of  arrivals. 

The  bottom  panel  of  Figure  7.1  shows  the  three  rays  predicted  by  IASP91 
to  exist  at  26°.  These  rays  form  a  classic  triplication  of  a  ray  turning  above 
the  discontinuity,  a  ray  reflecting  from  the  discontinuity,  and  a  ray  diving 
beneath  the  discontinuity.  In  this  case  the  diving  wave  arrives  first,  about 
three  seconds  ahead  of  and  1  s/deg  steeper  than  the  shallow  turning  ray.  The 
reflected  ray  arrives  just  0.3-0. 4  s  behind  the  shallow  turning  ray,  steeper  by 
about  0.3  s/deg.  So  while  it  is  unlikely  to  be  possible  to  discriminate  be¬ 
tween  the  two  shallower  waves  using  an  array  station,  it  may  be  possible  to 
discriminate  the  diving  ray  from  the  other  two. 

A  travel  time  curve  for  regional  distances  is  shown  in  Figure  7.2.  The 
curve  was  created  using  the  taup  methodology  of  Buland  &  Chapman  (1983) 
and  the  IASP91  velocity  model.  We  show  the  curve  in  reduced  space  in  which 
a  slowness  of  11  s/deg  was  used  to  offset  the  move-out.  The  cusps  associ¬ 
ated  with  the  transition  zone  phases  are  labeled  A  thru  F  and  are  ordered 
in  terms  of  decreasing  ray  parameter.  The  AB  branch  turns  above  the  410, 
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the  BC  branch  reflects  off  the  410,  the  CD  branch  turns  in  between  the  410 
and  the  660,  the  DE  branch  reflects  off  the  660,  and  the  EF  branch  turns 
below  the  660.  Points  E  and  C  are  predicted  to  have  large  amplitude  as  they 
represent  the  the  angle  of  critical  reflection  for  rays  incident  on  the  410  and 
660,  respectively.  Technically,  rays  reflected  at  precritical  angles  extend  the 
curves  at  points  C  and  E  back  towards  smaller  distances,  however  these  rays 
are  expected  to  have  small  amplitude.  Likewise,  diffraction  is  expected  to 
extend  the  curves  at  points  B  and  D  to  greater  distances,  but  again  ampli¬ 
tude  is  expected  to  drop  off  quickly,  especially  at  short  periods. 

Clearly,  different  velocity  models  lead  to  different  travel  time  curves.  For 
instance,  increasing  the  depth  to  the  410  km  discontinuity  delays  the  BC 
branch  in  time  and  shifts  it  to  larger  distances.  However,  for  reasonable 
velocity  models  the  general  character  of  the  curve  is  unlikely  to  change  fun¬ 
damentally,  and  the  IASP91  predictions  can  be  used  to  gain  insight  on  the 
observability  of  transition  zone  phases.  The  most  promising  distance  range 
to  distinguish  and  identify  triplicated  waves  in  observed  data  is  likely  at 
13°  —  15°.  Here  the  Pbc/cd  phases  are  expected  to  arrive  significantly  later 
and  less  steeply  than  Pab  and  shallower  phases,  and  to  have  relatively  large 
amplitude.  A  similar  situation  exists  at  distances  of  17°  —  18°  for  waves 
grazing  the  660  km  discontinuity.  The  third  promising  distance  is  around 
27°  —  28°  near  the  D  cusp  in  which  Pcd/de  rays  arrive  with  a  large  time 
separation  relative  to  Pef,  though  with  similar  ray  parameter. 

2.3  Overview  of  Previous  Work 

Many  previous  studies  have  used  seismic  arrays  to  study  upper  mantle  struc¬ 
ture.  When  the  first  teleseismic  arrays  were  deployed  in  the  early  1960’s,  it 
was  quickly  realized  that  the  direct  observations  of  dt/dA  (ray  parameter) 
provided  by  arrays  were  superior  to  those  obtained  with  the  more  common 
method  of  smoothing  and  differentiating  a  travel  time  vs.  distance  curve. 
The  derivative  information  ( dt/dA  as  a  function  of  A)  was  necessary  to  ob¬ 
tain  a  radial  velocity  model  via  the  Herglotz-Weichert  method.  One  of  the 
earliest  such  studies  was  that  of  Niazi  &  Anderson  (1965).  These  authors 
used  data  recorded  at  the  Tonto  Forest  array  in  Arizona  to  measure  dt/dA 
for  70  shallow  events  occurring  at  distances  of  10°  —  30°.  They  noticed  abrupt 
changes  in  dt/dA  at  distances  of  17°  and  24°  and  inferred  discontinuities  to 
exist  at  depths  of  320  km  and  640  km  respectively.  The  example  seismograms 
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presented  in  Niazi  &  Anderson  (1965)  show  clear  later  arrivals,  but  the  au¬ 
thors  were  unable  to  include  these  observations  in  the  modeling  because  of 
uncertainties  in  the  depths  of  the  sources.  Presumably,  the  authors  were 
worried  about  misinterpreting  depth  phases  as  triplicated  arrivals.  However, 
the  potential  value  of  these  arrivals  was  recognized,  with  the  final  sentence 
of  Niazi  &  Anderson  (1965)  stating,  “Clearly,  later  arrivals  must  be  used  to 
determine  details  of  transition  zones” . 

A  closely  related  study  was  published  two  years  later  using  a  similar 
data  set  recorded  at  the  Tonto  Forest  array  (Johnson,  1967).  In  this  case, 
the  author  was  able  to  observe  and  model  multiple  transition  zone  arrivals, 
stating  that  “A  general  conclusion  of  this  study  is  that  later  arrivals  pro¬ 
vide  crucial  data  for  eliminating  much  of  the  non-uniqueness  in  the  velocity 
structure”,  and  reporting  evidence  for  strong  gradients  (discontinuities)  at 
400  km  and  650  km.  The  primary  difference  in  the  two  studies  was  that 
while  Niazi  &  Anderson  (1965)  used  data  from  the  standard  10  km  aper¬ 
ture  cross  configuration  at  Tonto,  Johnson  (1967)  used  a  version  of  Tonto 
augmented  with  8  seismometers  that  yielded  an  effective  aperture  of  about 
300  km.  Most  likely,  the  increased  ray  parameter  resolution  provided  by  the 
increased  aperture  contributed  to  the  ability  of  Johnson  (1967)  to  confidently 
identify  and  model  the  triplicated  arrivals. 

A  thorough  study  of  the  upper  mantle  beneath  Australia  was  published 
in  1974  using  data  recorded  from  the  Warramunga  (WRA)  seismic  array 
(Simpson  et  ah,  1974).  The  authors  measured  the  ray  parameter  of  494 
P  wave  arrivals  from  sources  located  at  distances  of  10°  —  30°.  WRA  is  a 
medium-sized  array,  aperture  of  about  20  km,  and  the  authors  reported  a 
precision  of  better  than  0.1  s/deg  in  slowness  determination.  They  found 
pronounced  later  arrivals  at  21°  —  24.5°,  which  they  associated  with  a  dis¬ 
continuity  at  650  km,  and  large  amplitude  later  arrivals  at  12°  —  18°,  which 
they  associated  with  a  discontinuity  near  400  km.  Even  so,  they  commented 
on  the  non- uniqueness  of  the  problem  and  avoided  a  direct  inversion  of  their 
observations;  instead  they  used  forward  modeling  to  determine  a  consistent 
model  (SMAK  I)  and  properly  identify  the  later  arrivals. 

Similar  studies  involving  two  other  medium  aperture,  United  Kingdom 
style  arrays  were  carried  out  over  the  next  few  years.  Ram  &  Mereu  (1977) 
analyzed  data  from  over  350  earthquakes  at  distances  of  14°  —  36°  from  the 
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Gauribidanur  array  (GBA)  in  India.  The  authors  found  later  arrivals  and 
abrupt  changes  in  dt/dA  of  the  first  arrivals  that  were  consistent  with  dis¬ 
continuities  at  400  km  and  650  km.  They  also  discussed  the  possibility  of 
lateral  variation  in  MTZ  structure,  dividing  their  data  set  into  fourths  based 
on  backazimuth  and  inverting  for  separate  velocity  models.  Ram  et  al.  (1978) 
analyzed  data  from  over  100  earthquakes  recorded  at  the  Yellowknife  Array 
(YKA)  in  Canada  and  found  similar  results,  although  the  later  arrivals  were 
not  as  clearly  defined.  Transition  zone  multipathing  has  also  been  observed 
in  data  from  LASA  in  Montana  (see  Figure  6  of  Filson  (1975))  and  NORSAR 
in  Norway  (Ringdal,  1981),  although  we  did  not  find  any  dedicated  studies 
in  the  mainstream  geophysical  literature. 

In  more  recent  years,  the  trend  has  been  to  use  regional  seismic  networks 
(as  opposed  to  medium  aperture  arrays)  to  make  record  sections  of  MTZ  mul¬ 
tipathing  at  distances  of  10°  —  30°  for  a  relatively  small  number  (as  opposed 
to  a  large  number)  of  sources.  For  instance,  Walck  (1984)  used  data  from 
ten  earthquakes  recorded  in  Southern  California,  to  constrain  MTZ  struc¬ 
ture  beneath  the  Gulf  of  California;  Lefevre  &  Hclmberger  (1989)  modeled 
data  from  seven  earthquakes  recorded  at  seismometers  across  North  America 
to  infer  upper  mantle  structure  beneath  the  Canadian  Shield;  Ryberg  et  al. 
(1998)  used  data  from  approximately  a  dozen  peaceful  nuclear  explosions 
in  the  former  Soviet  Union  to  image  the  MTZ  in  northern  Eurasia;  Zhao 
et  al.  (1999)  used  data  from  eight  earthquakes  recorded  across  a  temporary 
seismometer  network  in  South  Africa;  Song  et  al.  (2004)  modeled  waveform 
data  generated  from  approximately  8  earthquakes  off  the  coast  of  Oregon 
and  Washington  and  recorded  at  regional  networks  throughout  the  Western 
U.S. 

Hence,  the  work  described  in  this  report  is  somewhat  of  a  “retro”  study, 
having  more  in  common  with  work  done  in  1960’s  and  1970’s  when  seismic 
arrays  were  relatively  new.  This  is  warranted  because  our  primary  interest  is 
not  in  inferring  Earth  structure,  but  in  understanding  wave  propagation  in 
South  Asia  in  the  context  of  verification  seismology.  The  primary  challenge 
is  that  the  array  considered  here  (CMAR)  was  designed  as  a  regional  ar¬ 
ray,  with  an  aperture  approximately  half  that  of  the  medium  aperture  arrays 
(WRA,  YKA,  GBA)  that  were  successfully  in  the  past  studies  of  upper  man¬ 
tle  multipathing  described  above.  Therefore,  we  can  expect  poorer  slowness 
resolution  and  potentially  more  difficulty  in  identifying  later  arriving  phases. 
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Array  Processing  Overview 

3.1  Description  of  Array  Processing  Techniques 

The  array  processing  software  used  in  this  study  was  previously  developed  by 
the  PI  under  the  acronym  GAP  (Generic  Array  Processing)  (Koper,  2001) 
and  has  been  slightly  modified  for  the  current  project.  The  core  problem  of 
slowness  inference  is  carried  in  the  time  domain  by  repeated  beam  formation 
over  a  2D  Cartesian  grid  of  slowness  vectors.  The  beam  giving  the  highest 
mean-square  power  over  the  prescribed  time  window  is  chosen  as  the  optimal 
beam,  and  the  corresponding  slowness  vector  is  chosen  as  the  optimal  slow¬ 
ness  vector.  This  technique  is  sometimes  referred  to  as  beampacking  and  is 
roughly  the  time-domain  equivalent  of  traditional  fk  techniques  (Schweitzer 
et  ah,  2002).  We  use  it  because  it  implicitly  averages  over  frequency,  while 
being  very  localized  in  time.  Therefore,  like  VESPA  analysis  (Davies  et  al., 
1971),  it  is  well-suited  for  analyzing  multiple  body  waves  that  arrive  in  short 
amount  of  time. 

The  user  is  required  to  enter  a  time  window,  corner  frequencies  and  num¬ 
ber  of  poles  for  a  bandpass  filter,  a  type  of  beam,  and  spacing  constraints 
on  the  slowness  grid.  In  this  study  we  commonly  use  3-polc  Butterworth 
bandpass  filters,  and  phase-stack  weighted  (PSW)  beams.  The  PSW  beams 
amplify  coherent  energy  but  induce  less  waveform  distortion  than  nth  root 
beams  (Schimmel  &  Paulssen,  1997).  We  find  that  this  leads  to  more  pre¬ 
cise  slowness  estimates.  The  traces  are  also  resampled  from  0.05  s  to  0.01  s 
to  increase  precision.  Each  array  element  is  weighted  equally  during  beam 
formation  and  no  special  allowance  is  made  for  the  local  noise  structure.  We 
experimented  with  inferring  3D  slowness  vectors  that  account  for  elevation 
differences  among  the  array  elements,  but  found  no  compelling  advantages. 
This  is  expected  because  CMAR  is  relatively  flat,  having  a  maximum  el¬ 
evation  difference  of  only  0.066  km,  compared  to  a  horizontal  aperture  of 

10.1  km  (see  Bokelmann  (1995)  for  a  discussion  of  topographic  effects  on 
slowness  inference). 

When  necessary,  error  estimates  are  made  using  a  bootstrap  technique 
(Tichelaar  &  Ruff,  1989)  that  generates  pseudo  arrays  by  resampling  the 
elements  with  replacement,  repeating  the  grid  search,  and  generating  a  pop¬ 
ulation  of  pseudo-solutions.  A  covariance  matrix  is  estimated  from  this  pop- 
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ulation  and  then  converted  into  a  95%  confidence  ellipse.  Also,  individual 
standard  errors  for  the  ray  parameter  and  backazimuth  are  obtained  from  the 
population  of  pseudo  solutions.  We  use  this  method,  as  opposed  to  comput¬ 
ing  some  sort  of  F-statistic,  because  our  definitions  of  power  generally  arise 
from  non-linear  beam  formation  techniques  that  are  not  necessarily  expected 
to  follow  standard  probability  distributions. 

We  also  operate  the  core  slowness  inference  process  in  a  sliding  window 
mode.  This  allows  for  the  generation  of  ray  parameter  and  backazimuth  time 
series,  in  addition  to  beam  power  and  coherence  as  a  function  of  time.  Movies 
that  show  power  as  a  function  of  ray  parameter,  backazimuth,  and  time  can 
also  be  output.  Our  measure  of  coherence  is  the  stack  of  the  instantaneous 
phases  of  the  aligned  traces,  which  varies  from  1  for  perfect  coherence  to  0 
for  perfect  incoherence  (Schimmcl  &  Paulssen,  1997).  In  this  work,  the  data 
are  generally  filtered  around  1.0  Hz,  and  we  have  found  a  window  length  of 
2  s  provides  good  results.  The  sliding  window  is  shifted  in  increments  of 
0.2  s,  giving  a  90%  overlap.  We  commonly  do  not  obtain  individual  error 
estimates  when  the  slowness  estimation  process  is  run  in  a  sliding  mode. 

3.2  Characteristics  of  CMAR 

3.2.1  Array  Geometry 

The  short-period  Chiang  Mai  array  consists  of  18  elements  arranged  in  a 
roughly  circular  geometry  (Figure  7.3).  The  resulting  153  inter-element  sep¬ 
arations  range  from  1.46  km  to  10.07  km,  reflecting  a  design  goal  of  detecting 
relatively  small  magnitude,  regional  distance  seismic  events.  Figure  7.3  also 
shows  the  CMAR  response  function  (as  defined  by  Rost  &  Thomas  (2002)) 
at  two  representative  frequencies.  Owing  to  the  circular  geometry  of  CMAR, 
the  ARF  is  very  symmetric  azimuthally,  and  can  be  well-approximated  as  a 
ID  function  of  ray  parameter.  At  1  Hz  the  half-width  of  the  center  peak  is 
about  6  s/deg,  while  at  3  Hz  the  half-width  is  about  2  s/deg.  In  both  cases, 
prominent  sidelobes  fall  outside  of  the  slowness  range  for  regional  distance 
P  waves,  and  spatial  aliasing  is  not  expected  to  be  a  problem. 
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3.2.2  Properties  of  the  Ambient  Noise  Field 

We  examined  the  ambient  noise  field  at  CMAR  using  data  recorded  during 
950  randomly  chosen,  4-sec  long  time  windows.  For  each  time  window  we 
calculated  a  2D  slowness  grid  with  the  horizontal  components  of  slowness 
varying  from  -40  s/deg  to  40  s/deg  in  increments  of  0.2  s/deg.  Beams  were 
computed  using  phase-stack  weighting  on  traces  bandpass  filtered  from  0.67- 
1.33  Hz  and  resampled  to  0.01  s.  Each  grid  was  then  normalized  by  its  peak 
amplitude  before  being  included  in  a  grid  stack.  Therefore,  our  technique 
does  not  recognize  the  absolute  level  of  noise  at  a  given  time,  but  only  its 
distribution  in  slowness  space. 

The  1-Hz,  time-averaged  ambient  noise  field  at  CMAR  possesses  three 
distinct,  coherent  peaks  that  are  robust  from  year-to-year  (Figures  7.4,  7.5, 
and  7.6).  The  most  prominent  occurs  to  the  southwest  at  a  backazimuth 
near  220  degrees  and  an  apparent  velocity  of  3. 5-4.0  krn/s;  the  next  largest 
occurs  at  apparent  velocities  above  15-16  krn/s,  generally  from  the  south 
although  backazimuth  is  poorly  determined  because  of  the  high  apparent  ve¬ 
locities;  the  smallest  peak  occurs  almost  due  East,  with  apparent  velocities  of 
3. 5-4.0  krn/s.  Importantly,  the  region  of  most  interest  in  a  monitoring  sense, 
basically  everything  arriving  from  the  North,  is  relatively  quiet  at  most  times. 

These  initial  observations  were  interesting  enough  that  we  did  a  compre¬ 
hensive  study  of  the  ambient  noise  at  CMAR.  The  full  analysis  is  described 
in  Koper  &  de  Foy  (2008),  which  is  attached  to  this  report  as  an  Appendix, 
and  here  we  just  mention  the  main  results: 

•  At  frequencies  above  1.4  Hz  the  noise  at  CMAR  is  isotropic  and  diffuse. 

•  At  frequencies  lower  than  1.4  Hz  the  noise  at  CMAR  is  strongly  par¬ 
titioned  by  apparent  velocity  into  two  categories:  teleseismic  P  wave 
energy  with  apparent  velocities  higher  than  25  km/s  (ray  parameters 
of  0. 0-5.0  s/deg)  and  higher  mode  Rayleigh  energy  with  apparent  ve¬ 
locities  near  4.0  km/s. 

•  The  Rayleigh  noise  is  further  partitioned  by  direction,  with  the  strongest 
signal  arriving  from  the  Bay  of  Bengal  at  backazimuths  of  180-255°.  A 
secondary  peak  in  the  Rayleigh  noise  occurs  in  the  direction  of  the 
South  China  Sea  at  backazimuths  of  80-120°. 
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•  The  Rayleigh  noise  is  strongly  seasonal  with  annual  variations  of  10- 
15  dB  in  power.  For  a  magnitude  scale  such  as  MS(Vmax)  (Russell, 
2006;  Bonner  et  al.,  2006)  that  depends  directly  on  the  logarithm  of 
seismic  amplitude,  our  observations  imply  that  the  detection  threshold 
at  CMAR  varies  by  0.5-0.75  magnitude  units  according  to  direction 
and  time  of  year. 

•  Like  the  Rayleigh  noise,  the  P  noise  observed  at  CMAR  is  seasonal. 
It  has  an  annual  power  variation  of  5-10  dB  and  there  are  several  geo¬ 
graphic  regions  that  could  act  as  sources:  the  western  Atlantic  Ocean 
near  the  coast  of  northern  Brazil  may  contribute  PKP  energy,  the  Pa¬ 
cific  Ocean  just  north  of  New  Guinea  may  contribute  PcP  energy,  and 
central  portions  of  the  North  Pacific  may  contribute  P  waves  that  turn 
in  the  lower  mantle. 

3.2.3  Signal  Coherence  vs.  Noise  Coherence 

We  also  examined  the  relative  coherence  of  signal  and  noise  as  a  function  of 
frequency  and  interstation  spacing,  using  data  from  950  events  that  occurred 
mainly  to  the  northwest  of  CMAR  (a  complete  description  of  these  events  is 
given  in  a  later  section).  We  used  the  following  procedure  for  each  event: 

1.  The  unfiltered  traces  were  resampled  to  0.01  s  and  aligned  on  the  first 
arrival  with  a  standard  cross-correlation  algorithm. 

2.  One  of  six  frequencies  bands  (centered  on  0.2  Hz,  0.5  Hz,  1.0  Hz,  2.0  Hz, 
3.0  Hz,  and  4  Hz)  was  chosen  and  a  narrow,  3-pole  Butterworth  band¬ 
pass  Liter  was  applied  to  each  trace. 

3.  Correlation  coefficients  for  the  153  distinct  element  pairs  were  com¬ 
puted  for  a  10  s  window  starting  2  s  before  the  first  arrival  (signal), 
and  a  10  s  window  starting  1  minute  before  the  first  arrival  (noise). 

The  approximately  145,000  correlation  coefficients  in  each  frequency  band 
(for  both  signal  and  noise  time  windows)  were  then  sorted  into  inter-element 
distance  bins  with  widths  of  1  km.  The  means  and  standard  deviations  of 
each  bin  are  shown  in  Figure  7.7  for  the  six  frequency  bands. 

Technically,  the  correlation  coefficients  should  be  Z-transformed  before 
the  mean  and  standard  deviations  are  computed  because  of  the  non-Gaussian 
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nature  of  the  correlation  coefficient  distribution  (van  Decar  &  Crosson,  1990); 
however,  this  correction  will  not  alter  the  fundamental  features  shown  in 
Figure  7.7.  We  find  that  the  frequency  band  giving  the  largest  gap  between 
signal  coherency  and  noise  coherency  is  centered  at  1  Hz,  and  therefore  this 
band  should  yield  the  highest  gain  during  beam-formation.  In  this  band,  the 
noise  is  almost  perfectly  uncorrelated  across  the  entire  array,  while  the  signal 
maintains  an  average  correlation  coefficient  near  0.55.  Using  equation  6  from 
Mykkcltveit  et  al.  (1983)  this  corresponds  to  an  average  gain  of  about  3.2, 
significantly  less  than  the  classic  gain  of  4.2  that  would  be  expected  if 
the  signal  were  perfectly  correlated.  It  is  important  to  note  that  the  preferred 
frequency  range  near  1  Hz  that  we  infer  is  potentially  only  valid  for  regional 
distance  events  to  the  northwest  of  CMAR.  As  shown  in  Figures  7.4,  7.5, 
and  7.6,  the  1  Hz  noise  in  other  geographical  quadrants  (and  at  small  ray 
parameters  in  general)  is  often  significantly  coherent.  Therefore,  we  might 
expect  the  largest  gap  between  signal  and  noise  coherence  to  be  at  a  higher 
frequency  band  for  these  directions. 
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4  Array  Analysis  of  CMAR  Data 

4.1  Chinese  Nuclear  Tests 


We  first  examined  data  from  three  Chinese  nuclear  explosions  that  were 
recorded  at  CMAR.  These  events  have  relatively  simple  source  time  functions 
and  no  depth  phases  that  could  interfere  with  later  arriving  direct  phases. 
High  quality  locations  for  these  three  events  are  given  by  Fisk  (2002)  as: 


Date 

Time 

Lat. 

Lon. 

rnh 

Distance 

Backazimuth 

1995/05/15 

04:05:59.38 

41.5545 

88.7516 

6.1 

24.65 

341.5 

1995/08/17 

00:59:59.35 

41.5402 

88.7533 

6.0 

24.64 

341.5 

1996/06/08 

02:55:59.37 

41.5780 

88.6875 

5.9 

24.69 

341.4 

where  distance  and  backazimuth  are  calculated  relative  to  the  beam  reference 
point.  At  this  distance  IASP91  predicts  three  arrivals  associated  with  the 
660  km  discontinuity.  The  first  is  an  Pef  wave  that  turns  beneath  the  660 
with  a  ray  parameter  of  9.1  s/deg;  about  1.3  s  later  a  Pcd  arrival  that  turns 
above  the  660  with  a  ray  parameter  of  10.3  s/deg  is  expected;  and  about 
1.1  s  later  a  Pee  wave  that  reflects  off  of  the  660  with  a  ray  parameter  of 
9.7  s/deg  is  expected. 

Record  sections  and  slant  stacks  for  these  three  events  are  shown  in  Fig¬ 
ures  7.8,  7.9,  and  7.10  respectively.  In  each  case  there  appears  to  be  three 
peaks  in  slowness  space  that  could  plausibly  be  associated  with  the  three  ex¬ 
pected  arrivals.  This  is  especially  true  for  the  May  15,  1995  explosion  shown 
in  Figure  7.8.  It  also  appears  that  the  first  arrivals  in  each  slantstack  have 
the  largest  ray  parameter,  perhaps  implying  that  for  this  path  the  depth  to 
the  660  km  is  slightly  larger  than  in  IASP91  and  so  the  Pcd  phase  arrives 
in  advance  of  Pef-  However,  each  peak  is  smeared  and  significantly  inter¬ 
feres  with  the  others.  Considering  that  these  three  events  are  essentially 
colocated,  it  is  disappointing  that  the  three  slantstacks  are  not  more  similar. 
Partially  this  is  due  to  this  specific  distance  being  near  a  branch  intersec¬ 
tion  on  the  traveltime  curve  (Figure  7.2),  however  it  is  also  indicative  of  the 
small  aperture  of  CMAR  as  discussed  in  the  introduction.  More  sophisticated 
techniques  (i.e.,  Capon,  1969;  Schmidt,  1986;  Shumway  et  ah,  2008)  could 
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potentially  provide  more  accurate  estimates  of  slownesses  for  triplicated  ar¬ 
rivals;  however,  all  slowness  estimators  are  likely  to  have  problems  when  the 
signals  arrive  from  similar  directions  at  similar  times  (Shumway  et  al.,  2008). 


4.2  Events  from  the  NEIC  Catalog  (1995-2004) 

Our  original  data  set  consisted  of  over  1,700  events  that  occurred  from  1995- 
2004  at  distances  of  13°  —  30°  from  CMAR.  We  selected  only  those  events 
that  had  primarily  a  continental  path  and  those  for  which  a  P  wave  arrival 
at  CMAR  was  reported  in  the  NEIC.  As  we  began  processing  these  data  we 
noticed  that  many  had  large  travel  time  residuals  relative  to  IASP91,  and  in 
some  cases  the  nominal  arrival  time  reported  in  the  NEIC  catalog  did  not 
match  a  noticeable  feature  in  the  data.  Therefore,  in  order  to  improve  the 
overall  quality  of  the  data  set  we  decided  to  cull  events  with  poor  quality 
locations. 

To  classify  which  events  to  keep  and  which  to  eliminate  we  consulted  the 
updated  EHB  catalog  of  re-processed  ISC  locations  (Engdahl  et  ah,  1998). 
We  accepted  those  events  from  our  original  data  set  for  which  (1)  the  isol 
field  was  HEQ,  DEQ,  or  FEQ,  (2)  there  was  no  x  in  the  iseq  field,  and  (3) 
the  ahyp  field  was  not  Z.  This  reduced  our  data  set  to  955  events,  but  led 
to  a  more  consistent  and  high  quality  data  set.  We  then  replaced  the  orig¬ 
inal  NEIC  hypocenters  with  the  EHB  hypocenters.  The  distribution  and 
statistics  of  this  data  set  are  shown  in  Figure  7.11.  It  is  complete  down  to  a 
magnitude  of  about  4.7  uif,.  Most  of  the  seismicity  forms  an  arcuate  belt  rep¬ 
resentative  of  the  ongoing  collision  between  India  and  Eurasia.  Also,  most  of 
the  seismicity  is  shallow,  with  the  main  exception  being  intermediate  depth 
events  in  the  Hindu-Kush  zone  in  Afghanistan.  An  example  record  section 
of  array  beams  is  shown  in  Figure  7.12. 


4.2.1  Travel  Time  Residuals  of  First  Arrivals 

For  each  event  we  visually  inspected  the  traces  and  eliminated  those  that  had 
glitches,  null  segments,  or  electronic  noise.  We  then  filtered  the  data  using 
a  3-pole  Butterworth  bandpass  with  corners  of  0.67  Hz  and  1.33  Hz,  formed 
linear  beams  using  the  expected  slowness  vector,  and  hand  picked  the  arrival 
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times  of  the  first  arriving  waves.  This  was  partially  done  for  quality  control, 
but  also  because  these  times  were  needed  to  align  beams  prior  to  stacking 
them  in  bins  of  distance  or  azimuth.  Our  initial  processing  indicated  that 
theoretical  times,  even  for  EHB  locations,  were  not  accurate  enough  to  pro¬ 
duce  coherent  beam  stacks. 

The  corresponding  travel  time  residuals  are  shown  geographically,  and 
averaged  in  bins  of  distance,  in  Figure  7.13.  The  distribution  is  nearly  Gaus¬ 
sian  with  a  mean  of  -0.11  s  and  a  standard  deviation  of  1.64  s.  We  attribute 
the  remarkably  small  bias  to  the  high-quality  of  the  EHB  locations  and  the 
appropriateness  of  IASP91  as  an  average  velocity  model  for  the  region.  Co¬ 
herent  geographical  patterns  appear  in  the  residuals,  such  as  the  patch  of 
near-zero  values  for  events  in  Pakistan  and  the  strongly  negative  signal  for 
events  in  Tibet;  however,  it’s  unclear  if  these  patterns  are  related  to  true 
3D  variations  in  geology  or  simply  reflect  an  inaccuracy  in  the  ID  reference 
model  (note  the  distance  dependence  to  the  travel  time  residuals).  Obvi¬ 
ously,  data  from  many  more  stations  and  earthquakes  would  be  required  to 
generate  the  crossing  paths  required  for  a  tomographic  study. 


4.2.2  Slowness  Residuals  of  First  Arrivals 

A  previous  study  of  slowness  anomalies  at  CMAR  caused  by  near  receiver 
heterogeneities  found  a  median  ray  parameter  residual  of  —1.35  ±0.58  s/deg 
and  a  median  backazimuth  residual  of  7.3  ±  12.4°  (Bondar  et  ah,  1999).  This 
level  of  inaccuracy  is  large  enough  to  cause  confusion  among  various  upper 
mantle  branches  of  the  travel  time  curve,  especially  near  cusps  (Figure  7.2). 
However,  those  authors  also  found  that  the  slowness  anomalies  depended  sig¬ 
nificantly  on  incoming  direction  and  the  values  reported  above  are  an  average 
over  all  regions  of  slowness  space.  Therefore,  as  a  further  check  on  data  qual¬ 
ity,  but  more  importantly  to  assess  the  accuracy  of  CMAR-based  slowness 
estimates  for  regional  distance  events  from  the  northwest,  we  estimated  2D 
slowness  vectors  for  all  of  the  first  arrivals  in  our  data  set. 

We  filtered  the  data  with  a  3-pole  Butterworth  bandpass,  using  corner 
frequencies  of  0.67  Hz  and  1.33  Hz,  and  determined  the  optimal  slowness 
vector  for  each  event  by  doing  a  grid  search  over  beam  power  in  Cartesian 
slowness  coordinates.  The  beams  were  calculated  using  third  order  phase 
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stack  weighting  and  power  was  measured  over  a  4  s  window  that  started 
1  s  before  the  picked  first  arrival  time.  In  Figure  7.14  we  show  the  re¬ 
sulting  residuals  assuming  EHB  locations  and  the  IASP91  velocity  model. 
We  find  mean  residuals  of  —0.84  ±  1.14  s/deg  and  —0.28  ±  1.16  s/deg  for 
East  /West  and  North/South  slowness  respectively.  The  trend  is  remarkably 
constant  for  various  incoming  directions  and  therefore  is  probably  created 
by  some  sort  of  geological  heterogeneity  near  the  array.  In  general  though, 
the  ray  parameter  estimates  are  only  mildly  biased,  with  a  mean  residual  of 
0.18  ±  1.26  s/deg.  The  observed  ray  parameters  nicely  track  the  expected 
change  of  about  5  s/deg  as  distance  increases  from  13°  to  30°.  This  implies 
that  measured  ray  parameters  for  later  arrivals,  at  least  near  cusps  C  and  E, 
should  be  distinguishable  from  the  first  arrival  ray  parameters. 


4.2.3  Properties  of  Later- Arriving  Coda  Waves 

We  experimented  with  several  approaches  in  an  effort  to  make  images  of  co¬ 
herent  later  arriving  phases.  For  instance,  we  tried  stacking  vespagrams  that 
had  been  formed  at  theoretical  backazimuths  into  bins  according  to  distance, 
and  we  tried  coherent  and  incoherent  stacking  of  array  beams  for  specific 
source  regions.  In  general  though,  these  efforts  were  unsuccessful  and  we 
ultimately  settled  on  a  technique  that  simply  identified  significant  later  ar¬ 
rivals  on  an  event-by-event  basis.  We  then  binned  these  arrivals  by  distance 
and  looked  for  trends  in  slowness,  timing,  and  relative  amplitude. 

Our  technique  for  identifying  significant  later  arrivals  is  described  as  fol¬ 
lows.  For  each  event  we  performed  a  sliding  window  slowness  analysis,  start¬ 
ing  10  s  prior  to  the  picked  first  arrival  time  and  continuing  for  40  s.  The 
windows  were  two  seconds  long  and  were  shifted  in  increments  of  0.2  s,  giving 
a  90%  overlap.  Optimal  slowness  vectors  were  determined  for  each  window 
using  a  grid-search  method  similar  to  what  was  described  in  the  previous 
section  and  with  the  same  frequency  band.  For  each  event,  this  resulted  in 
four  time  series  with  equal  lengths  of  201  points:  beam  power,  coherence,  ray 
parameter,  and  backazimuth.  Example  output  for  this  procedure  is  shown 
in  Figure  7.15. 

We  then  tested  several  criteria  for  identifying  significant  arrivals  from 
these  time  series.  The  following  requirements  led  to  a  nice  balance  between 
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not  generating  too  many  spurious  arrivals  while  not  missing  legitimate  but 
small-amplitude  arrivals.  First,  we  multiplied  the  linear  beam  power  by  the 
square  of  the  coherence  and  applied  two  passes  of  a  nearest-neighbor  smooth¬ 
ing  filter.  Any  local  maximum  on  this  curve  was  then  accepted  as  an  arrival 
if  it  had  (1)  a  signal-to- noise  ratio  of  at  least  5  on  the  original  linear  beam, 
(2)  a  coherence  of  at  least  0.8,  and  (3)  a  ray  parameter  estimate  less  than 
16  s/deg.  For  the  event  shown  in  Figure  7.15,  which  is  typical,  this  proce¬ 
dure  resulted  in  seven  distinct  arrivals  as  shown  by  the  red  stars.  Modest 
changes  in  our  acceptance  criteria  led  to  small  or  no  changes  in  the  number 
of  significant  arrivals. 

Results  from  this  procedure  are  shown  in  Figure  7.16  for  a  group  of  shal¬ 
low  events  (depth  less  than  40  km)  that  occurred  at  distances  of  14.5°  — 15.5°. 
The  amplitude  of  each  later  arrival  is  scaled  according  to  the  event-specific 
maximum  and  then  plotted  as  a  function  of  ray  parameter,  backazimuth, 
and  lapse  time.  There  is  some  scatter  in  the  backazimuth  estimates  but 
they  clearly  cluster  around  the  value  for  the  first  arrival.  This  is  true  even 
for  the  very  small  amplitude  arrivals  that  are  only  10%-20%  as  large  as  the 
maximum.  We  also  found  this  to  be  true  for  the  other  distance  ranges  in 
our  data  set.  This  implies  that  most  of  the  P-coda  we  observe  is  created 
by  near-source  heterogeneity  and  multipathing.  While  some  near-receiver 
P  —  Rg  scattering  may  occur,  there  appears  to  be  very  little  near-receiver 
P  —  P  scattering. 

From  the  stack  of  normalized  beam  power  time  series  shown  in  Fig¬ 
ure  7.16,  it  is  clear  that  there  are  two  distinct  pulses  of  energy  for  events 
occurring  at  distances  around  15°.  As  shown  in  the  upper  right  panel,  the 
secondary  energy  arrives  about  6-8  s  after  the  first  arrival  and  has  ray  param¬ 
eters  2-3  s/deg  smaller  than  the  first  arrival.  These  properties  are  consistent 
with  Pbc  arrivals  that  have  reflected  off  of  the  410-km  discontinuity  and/or 
Pcd  arrivals  that  have  turned  beneath  the  410-km  discontinuity.  It  is  the 
drop  in  ray  parameter  in  particular  that  distinguishes  these  arrivals  from 
pPab  depth  phases.  As  a  numerical  example,  consider  a  40  km  deep  event 
occurring  at  15°  from  CMAR  (such  deep  crustal  events  are  relatively  common 
in  Tibet).  The  following  arrivals  are  predicted  by  IASP91: 
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Phase 

Travel  Time  (s) 

Ray  Parameter  (s/deg) 

P  (AB) 

- 

13.1 

P  (asth.) 

0.20 

13.6 

P  (asth.) 

0.35 

13.5 

P  (CD) 

4.87 

11.1 

P  (BC) 

4.97 

11.2 

pP  (AB) 

8.00 

13.6 

pP  (asth.) 

8.26 

13.3 

pP  (asth.) 

8.28 

13.5 

pP  (CD) 

14.72 

11.1 

pP  (BC) 

14.74 

11.2 

So,  while  depth  phases  from  shallow  P  waves  could  plausibly  explain  the 
observed  differential  times,  they  cannot  explain  the  observed  drop  in  ray  pa¬ 
rameter.  Furthermore,  synthetic  seismograms  show  that  interfering  Pbc/cd 
waves  near  the  C  cusp  often  have  bigger  amplitudes  than  the  first- arriving 
Pab  wave.  Finally,  we  point  out  that  these  strong  later  arrivals  show  up 
at  distances  out  to  18°,  but  at  progressively  smaller  differential  times  (Fig¬ 
ure  7.17).  This  is  the  expected  behavior  for  Pcd  arrivals,  but  not  for  pPab 
depth  phases. 

The  mean  differential  time  between  the  largest  arrival  and  first  arrival  for 
the  15°  distance  bin  is  6.9  +/-  1.7  s  (after  removing  4  outliers).  This  is  about 
2  s  larger  than  predictions  from  IASP91  for  a  surface  focus,  implying  either 
that  Pab  is  too  fast,  Pbc/cd  is  too  slow,  or  some  combination  of  the  two. 
We  previously  found  evidence  that  Pab  waves  (the  first  arrivals)  are  fast  for 
this  area  (Figure  7.13);  in  fact,  the  mean  absolute  Pab  residual  for  these 
events  is  about  -2  s.  So  the  simplest  explanation  is  that  the  uppermost  man¬ 
tle  beneath  Burma  is  seismically  fast,  perhaps  due  to  extra  thick  lithosphere 
or  lid,  while  the  region  around  the  410  km  is  nicely  described  by  IASP91. 
Alternatively,  some  of  the  anomaly  in  Pab-Pbc/cd  differential  times  could 
be  caused  by  a  slightly  depressed  410  km  discontinuity  in  the  region.  A  de¬ 
pression  of  10  km  would  add  about  0.5  s  to  the  differential  time.  Given  this 
level  of  non-uniqueness,  formal  modeling  of  Earth  structure  is  not  warranted. 

In  Figure  7.18  we  show  depictions  of  significant  later  arrivals  for  distances 
of  19°  —  24°.  The  Pde/ef  branch  of  arrivals  is  expected  to  appear  5-6  s  af¬ 
ter  the  first  arrival  with  significantly  reduced  ray  parameter,  but  there  is  no 
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compelling  trend  in  the  observations.  This  may  be  related  to  the  sparseness 
of  data  at  distances  of  22°  —  24°  and  is  not  necessarily  an  indicator  of  non¬ 
standard  structure  near  the  660  km  discontinuity.  At  the  largest  distances  of 
25°  —  30°  (Figure  7.19)  we  again  see  evidence  of  multipathing  in  the  upper 
mantle.  Especially  at  distances  of  27°  —  29°  secondary  arrivals  appear  2-6  s 
after  the  first  arrival  with  significantly  increased  ray  parameter,  as  predicted 
for  Poe  arrivals.  However,  there  is  wide  variation  in  the  differential  times 
and  ray  parameters,  and  no  strong  constraints  can  be  placed  on  Earth  struc¬ 
ture. 

ft  is  worth  noting  that  at  the  largest  distances  of  28°  —  30°  the  stacked 
beam  power  envelopes  flatten  out  considerably,  with  many  strong  arrivals  oc¬ 
cupying  the  entire  time  window  of  30  s  after  the  first  arrivals  (Figure  7.19). 
Only  shallow  events  (depth  <  40  km)  are  included  in  these  stacks,  and  so  it 
cannot  be  that  the  arrivals  are  depth  phases  from  intermediate  depth  earth¬ 
quakes  in  the  Hindu-Kush  seismogenic  zone.  Instead,  the  nearly  flat  coda 
envelopes  are  likely  indicative  of  increased  near  source  scattering.  A  plausi¬ 
ble  (though  non-unique)  explanation  of  this  energy  is  scattering  of  Rg  to  P 
from  topographic  features  in  the  source  region.  The  relatively  slow  velocity 
of  Rg  provides  a  mechanism  for  generating  significant  P  energy  far  after  the 
first  arrival  with  only  a  single  scattering  event.  Evidence  exists  for  such  a 
mechanism  based  on  nuclear  explosion  data  from  Nevada  Test  Site  (Stead  & 
Helmberger,  1988;  Gupta  et  ah,  1991).  However,  for  this  mechanism  to  be 
viable  the  earthquakes  in  this  area  would  have  to  be  especially  shallow,  as 
Rg  excitation  drops  off  significantly  with  increasing  source  depth. 
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5  Conclusions  and  Recommendations 


Using  the  medium-aperture  array  in  Chiang  Mai,  Thailand  we  studied  the 
properties  of  P  coda  waves  for  a  group  of  955  events  that  occurred  at  dis¬ 
tances  of  13°  —  30°.  We  found  that  nearly  all  coherent  energy  in  the  first 
30  s  after  the  P  wave  arrives  along  the  great  circle  path.  This  indicates 
that  near  receiver  scattering,  which  would  generally  lead  to  large  slowness 
anomalies,  contributes  very  little  to  the  coda  waves.  Nearly  all  coherent  en¬ 
ergy  also  arrives  with  ray  parameters  within  2-3  s/deg  of  the  first  arrival.  The 
majority  of  this  energy  seems  to  be  associated  with  near-source  scattering, 
perhaps  Rg  into  P,  though  this  would  require  especially  shallow  earthquakes. 

Although  we  were  unable  to  produce  to  images  of  transition  zone  multi- 
pathing  with  various  beam-stacking  techniques,  we  did  observe  clear  evidence 
for  upper  mantle  multipathing  from  events  occurring  in  Tibet  at  distances 
of  14°  —  17°  from  CMAR.  We  observed  secondary  arrivals  associated  with 
ray  paths  that  either  dive  beneath  the  410  km  discontinuity  ( Pcd )  or  reflect 
off  of  it  (PBC),  that  had  anomalously  large  differential  times  relative  to  the 
first  arriving  Par  phase.  The  differential  times  were  too  large  by  about  2  s 
which  is  nearly  equal  to  the  absolute  Par  residuals  determined  from  EHB 
locations.  Therefore,  while  a  depressed  410  km  discontinuity  may  contribute 
to  the  anomalous  differential  travel  times,  it  is  more  likely  that  the  shallow 
mantle  beneath  Burma  is  anomalously  fast,  perhaps  because  of  a  thickened 
lid. 


In  general,  the  ray  parameter  and  travel  time  observations  made  in  this 
study  were  useful  in  identifying  the  nature  of  the  regional  distance  P  coda 
waves,  but  were  not  of  high  enough  quality  to  place  meaningful  constraints  on 
Earth  structure.  A  larger  aperture  array  (i.e.,  YKA  style  of  approximately 
20  km)  would  improve  the  slowness  resolution  and  perhaps  lead  to  better 
results,  however  the  existence  of  near-source  scattered  energy  would  still  be 
problematic.  In  terms  of  constraining  Earth  structure  from  upper  mantle 
triplications,  a  more  fruitful  approach  is  to  use  regional  networks  as  wide- 
aperture  arrays  for  a  relatively  small  number  of  well-located,  high-quality 
events.  This  gives  record  sections  across  large  segments  of  distance  with 
equivalent  source  signatures. 

The  most  interesting  work  accomplished  in  this  contract  involved  the 
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study  of  the  ambient  noise  field  at  CMAR,  in  which  we  found  strong  anisotropy 
in  both  apparent  velocity  and  backazimuth  at  frequencies  smaller  than  about 
1.4  Hz.  Noise  was  peaked  at  velocities  corresponding  to  teleseismic  P  waves 
and  higher  mode  Rayleigh  waves,  while  the  ring  of  slowness  space  corre¬ 
sponding  to  regional  distance  P  waves  was  relatively  quiet.  The  Rayleigh 
noise  arrived  preferentially  from  the  southwest  and  the  east,  while  northern 
directions  corresponding  to  the  main  Asian  landmass  were  quiet.  All  the 
noise  sources  were  well  correlated  with  ocean  wave  heights,  showing  strong 
seasonal  variation  in  power.  Hence,  detection  thresholds  at  CMAR  for  tele- 
seismic  body  waves  and  regional  surface  waves  are  seasonally  and  direction- 
ally  dependent. 
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7  Figures 

7.1  Transition  Zone  Ray  Paths 


Ray  paths  through  IASP91  for  (top)  ray  parameters  from  8.0-15.0  s/deg  in 
increments  of  0.05  s/deg,  (middle)  the  five  geometric  rays  predicted  to  exist 
at  17°,  and  (bottom)  the  three  geometric  arrivals  predicted  to  exist  at  26°. 
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7.2  Transition  Zone  Travel  Time  Curve 


IASP91  travel  time  curve  for  a  surface  focus,  reduced  by  11  s/deg.  Two  as- 
thenospheric  branches  merge  at  a  cusp  near  18.5°  and  we  define  this  energy 
as  Pn.  In  real  data  this  sort  of  energy  will  be  observable  to  larger  distances 
and  at  slightly  larger  ray  parameters  as  underside,  whispering  gallery  reflec¬ 
tions  from  the  Moho  contribute.  The  transition  zone  branches  are  shown  in 
red  and  labeled  A-F  in  order  of  decreasing  ray  parameter. 
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North/South  Slowness  (s/deg) 


7.3  Array  Response  of  CMAR 
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-40  -30  -20  -10  0 

Normalized  Power  [dB] 


(upper  left)  Station  geometry  of  CMAR  and  corresponding  array  response 
at  (lower  left)  1  Hz  and  (lower  right)  3  Hz.  (upper  right)  The  az- 
imuthally  averaged  array  response  at  1  Hz  and  3  Hz. 
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7.4 


Yearly  Averages  of  Ambient  Noise  at  CMAR 


1996 


2003 


1998 


2002 


1999 


997 


2001 


2000 


CMAR  noise  at  1  Hz  averaged  over  nine  years  (center)  and  averaged  in  one- 
year  bins.  Amplitudes  across  slowness  space  are  computed  using  phase-stack 
weighting  and  then  normalized  before  being  time-averaged. 
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7.5  Monthly  Averages  of  Ambient  Noise  at  CMAR 


November 


December 


January 


Winter 


October 


Fall 


Summer 


Spring 


March 


April 


CMAR  noise  at  1  Hz  from  a  nine-year  period  binned  according  to  month. 
Amplitudes  across  slowness  space  are  computed  using  phase-stack  weighting 
and  then  normalized  before  being  time-averaged. 
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7.6  Hourly  Averages  of  Ambient  Noise  at  CMAR 


18-19 


Local  Night 
(9  pm  -  9  am) 


Local  Day 
(9  am  -  9  pm) 


06-07 


12-13 


CMAR  noise  at  1  Hz  from  a  nine- year  period  binned  according  to  GMT  hour. 
Amplitudes  across  slowness  space  are  computed  using  phase-stack  weighting 
and  then  normalized  before  being  time-averaged. 
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7.7  Signal  and  Noise  Coherence  at  CMAR 


Sensor  Separation  (km)  Sensor  Separation  (km) 

Correlation  of  signal  (black)  and  noise  (red)  as  a  function  of  inter-element 
separation  for  six  frequency  bands.  The  signal  was  defined  using  a  10-s 
window  around  the  first  arrival  for  approximately  950  events  to  the  northwest 
of  CMAR,  while  the  noise  was  defined  using  a  10-s  window  centered  about 
a  minute  before  the  corresponding  first  arrival. 
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7.8  Chinese  Nuclear  Explosion  of  May  15,  1995 
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Record  section  and  slantstack  at  CMAR  for  the  Chinese  nuclear  explosion  of 
May  15,  1995.  The  data  were  filtered  with  a  3-polc  bandpass  at  0.67-2.67  Hz 
and  phase  weighted  stacking  of  order  3  was  used. 
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7.9  Chinese  Nuclear  Explosion  of  August  17,  1995 
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Record  section  and  slantstack  at  CMAR  for  the  Chinese  nuclear  explosion 
of  August  17,  1995.  The  data  were  filtered  with  a  3-pole  bandpass  at  0.67- 
2.67  Hz  and  phase  weighted  stacking  of  order  3  was  used. 
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7.10  Chinese  Nuclear  Explosion  of  June  8,  1996 
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Record  section  and  slantstack  at  CMAR  for  the  Chinese  nuclear  explosion  of 
June  8,  1996.  The  data  were  filtered  with  a  3-pole  bandpass  at  0.67-2.67  Hz 
and  phase  weighted  stacking  of  order  3  was  used. 
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7.11 


Locations  of  EHB  Events 


0 


34  70 


150 

Depth  (km) 


300 


Locations  of  the  955  events  from  the  EHB  catalog  that  we  considered  in  this 
study.  Circle  size  is  proportional  to  magnitude  and,  as  shown  in  the  inset, 
the  catalog  is  complete  down  to  a  magnitude  of  about  4.7  mb. 
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7.12  Example  Beam  Record  Section 


200  210  220  230  240  250 

Time  (s) 


260  270 


280 


Example  record  section  from  a  corridor  of  backazimuths  at  305°  —  310°.  The 
traces  are  phase-stack  weighted  beams  formed  at  the  expected  slowness,  and 
the  red  curve  is  the  predicted  travel  time  curve  for  IASP91.  The  times  have 
been  reduced  by  8  s/deg. 
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7.13  Travel  Time  Residuals  of  First  Arrivals 


-4  -3  -2-10  1 

Travel  Time  Residual  (s) 


Travel  time  residuals,  relative  to  IASP91,  for  the  EHB  data  set  considered  in 
this  study.  The  arrival  times  were  handpicked  and  used  to  facilitate  alignment 
of  array  beams. 
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7.14  Slowness  Residuals  of  First  Arrivals 


Sx  Residual  (s/deg) 


Sy  Residual  (s/deg) 


Distance  (deg)  Sx  (sec/deg) 


Slowness  analysis  of  first  arriving  phases  for  the  EHB  data  set.  The  two 
upper  panels  show  histograms  of  slowness  residuals  in  Cartesian  coordinates, 
the  lower  left  panel  shows  binned  ray  parameter  estimates  and  predictions, 
and  the  lower  right  panel  shown  full  vector  slowness  residuals  in  the  manner 
of  Bondar  et  al.  (1999). 
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7.15  Example  Sliding  Window  Slowness  Analysis 


Time  (s) 


Example  output  from  the  sliding  window  slowness  analysis  of  an  event  located 
15.1°  from  CMAR.  The  top  panel  shows  a  simple  linear  beam  formed  at  the 
expected  slowness  while  the  next  four  panels  show  beam  power,  coherence, 
ray  parameters,  and  backazimuth  as  a  function  of  time.  The  stars  indicate 
arrivals  picked  by  the  automated  algorithm  described  in  the  text. 
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Normalized  Power 


7.16  Properties  of  P  coda  at  Distances  Near  15° 
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1.0 


Coda  properties  for  83  shallow  events  that  occurred  at  distances  of  14.5°  — 
15.5°  from  CMAR.  The  upper  left  panel  shows  the  event  locations  while  the 
lower  left  panel  shows  a  beam  power  stack.  The  grey  area  represents  ±  one 
standard  deviation.  The  panels  on  the  right  show  the  ray  parameters  and 
backazimuths  for  significant  arrivals  within  the  coda,  as  determined  by  the 
algorithm  described  in  the  text. 
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7.17  P  coda  at  Distances  of  13°  —  18° 
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Panels  on  the  left  show  stacks  of  beam  power  from  events  in  various  distance 
bins.  The  grey  area  indicates  ±  one  standard  deviation.  Panels  on  the  right 
show  the  estimated  slowness  for  significant  arrivals  in  the  coda. 
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7.18  P  coda  at  Distances  of  19°  —  24° 
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Panels  on  the  left  show  stacks  of  beam  power  from  events  in  various  distance 
bins.  The  grey  area  indicates  ±  one  standard  deviation.  Panels  on  the  right 
show  the  estimated  slowness  for  significant  arrivals  in  the  coda. 
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Normalized  Power  Normalized  Power  Normalized  Power  Normalized  Power 


7.19  P  coda  at  Distances  of  25°  —  30° 
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Panels  on  the  left  show  stacks  of  beam  power  from  events  in  various  distance 
bins.  The  grey  area  indicates  ±  one  standard  deviation.  Panels  on  the  right 
show  the  estimated  slowness  for  significant  arrivals  in  the  coda. 
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9  Appendix 


9.1  Preprint  of  Koper  &  de  Foy  (2008) 

In  this  section  we  include  a  preprint  of  a  manuscript  that  arose  from  this 
contract: 

Koper,  K.D.,  and  B.  de  Foy,  2008.  Seismic  noise  from  Earth’s  deep  interior, 
BSSA ,  under  review. 

It  was  submitted  to  BSSA  on  May  19,  2008  and  is  currently  under  review. 
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i  Abstract 

i  We  present  an  analysis  of  seismic  noise  recorded  during  1995-2004  by  a  medium- 
s  aperture,  short-period  seismic  array  located  in  Chiang  Mai,  Thailand  (CMAR).  We 
«  calculated  frequency-wavenumber  spectra  for  nearly  1000  randomly  selected  time 
s  windows,  each  with  a  length  of  160  s.  At  frequencies  above  about  1.4  Hz  the  noise 
t  is  unorganized  and  wavenumber  spectra  are  isotropic  and  diffuse;  however,  at  lower 

i  frequencies  three  robust  wavenumber  peaks  exist.  Two  of  the  peaks  have  phase 
e  velocities  centered  near  4.0  km/s,  consistent  with  higher-mode  Rayleigh  waves, 
»  while  the  third  peak  has  much  higher  apparent  velocity  (>25  km/s),  consistent 
is  with  waves  that  have  interacted  with  Earth's  core  (PKP,  PcP).  All  three  peaks  are 

ii  strongly  seasonal  with  annual  power  variations  of  10-20  dB,  and  all  show  excellent 
«  correlation  in  their  putative  source  regions  with  ocean  wave  heights  derived  from 
is  TOP  EX/ POSEIDON  satellite  tracks.  To  the  best  of  our  knowledge,  this  is  the  first 
i«  time  such  a  high-velocity  component  of  seismic  noise  has  been  consistently  ohserved. 
is  The  presence  of  this  high-velocity  peak  raises  the  possibility  of  using  ambient  noise 
is  to  image  Earth  s  lower  mantle  and  core. 
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1  Introduction 


i*  Seismometer  arrays  have  been  used  to  study  seismic  noise  at  frequencies  span- 
x  ning  more  than  three  orders  of  magnitude,  from  as  low  as  0.002  Hz  (Rhie  & 
a  Romanowiez,  2004)  to  above  2  Hz  (Okada,  2003).  Arrays  are  exceptionally 
a  useful  in  noise  studies  because  they  provide  directional  constraints  on  prop- 
25  agating  energy  (Rest  &  Thomas,  2002).  Measurements  of  apparent  velocity' 
a  allow'  body  w'ave  components  of  noise  to  be  identified  and  distinguished  from 
a  surface  w'ave  components  (Toksoz  &  Lacces,  1968;  Haubrich  Sc  McCamy,  1969; 
2»  Lacces  et  al.,  1969)  and  allow  surface  wave  dispersion  curves  to  be  observed 
v  and  inverted  for  geologic  models  of  the  subsurface  (e.g.  Louie,  2001).  Mea- 
:»  surements  of  backazimuth  allow  specific  sources  of  seismic  noise  to  be  located 
»  (Bungum  et  al.,  1971;  Cessaro,  1994;  Schulte-Pelkum  et  al.,  2004)  and  even 
»  tracked  as  a  function  of  time  (Gerstoft  et  al.,  2006). 

jj  Here  we  present  an  analysis  of  seismic  noise  recorded  over  a  ten-year  period 
»  at  a  seismic  array'  in  Chiang  Mai,  Thailand.  Known  as  CMAR,  the  array  is  a 
»  primary  station  of  the  International  Monitoring  System  (IMS)  and  is  located 
m  at  far-regional  distances  from  nuclear  test  sites  in  China,  India,  Pakistan,  and 

*  North  Korea.  Part  of  the  motivation  for  this  study  is  to  determine  noise  char- 
»  acteristics  at  CMAR  that  might  affect  its  nuclear  monitoring  capability.  For 

*  instance,  detection  thresholds  at  CMAR  likely  depend  significantly  on  the  sea- 
se  son,  direction,  and  apparent  velocity.  More  generally,  the  noise  characteristics 
»  of  CMAR  are  also  interesting  because  they  illuminate  fundamental  dynamic 
4*  interactions  among  the  oceans,  atmosphere,  and  solid  Earth  (e.g.  Webb,  1998). 
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«  2  Characteristics  of  CMAR 

*i  CMAR  consists  of  18  vertical-component  seismometers  arranged  in  a  roughly 
45  circular  geometry  (Figure  la).  The  153  inter-element  separations  range  from 
44  1.46  km  to  10.07  km.  and  reflect  a  design  goal  of  detecting  teleseismic  and 
4*  far -regional  P  waves  at  frequencies  around  1  Hz.  The  corresponding  array  re- 
44  sponse  function  (Figure  lb)  is  azimuthally  symmetric  and  has  sidelobes  that 
47  are  small  and  well-separated  from  the  main-lobe.  Although  CMAR  is  con- 

44  sidered  a  medium-apert.ure  anay  by  the  seismic  monitoring  community,  the 
4*  inter-element  separations  are  small  enough  that  near  receiver  heterogeneities 
»  can  create  travel  time  anomalies  that  are  a  significant  fraction  of  the  move-out 
m  across  the  array.  Potentially,  this  can  lead  to  biased  estimates  of  the  slowness 
m  vector.  Examples  of  such  heterogeneities  observed  at  other  arrays  include  a 

53  dipping  Moho  (Tibuleac  &  Herrin,  1997),  intra- array  topography  ( Bokelmann. 

54  1995),  and  general  variations  in  geology  (Engdahl  k  Felix,  1971). 

54  Fbrtunately,  such  slowness  biases  have  been  found  to  be  relatively  mild  at 
m  CMAR  (Bondar  et  al.,  1999).  In  that  study,  the  authors  compared  observed 
57  and  predicted  slowness  vectors  at  CMAR  for  a  large  number  of  well-lc>cated 
54  earthquakes  and  found  a  median  raj'  parameter  anomaly  of  - 1 .35  ±  0.58  s/deg 
•»  and  a  median  backazimuth  anomaly  of  7.3  ±  12.43°.  Most  of  this  bias  can  be 
w  counteracted  by  applying  corrections  of  -1.622  s/deg  and  0.208  s/deg  to  ob- 
4i  servations  c>f  east-west  and  north-south  slowness  respectively.  The  bias  can  be 
4i  completely  eliminated  if  additional,  directionally  dependent,  corrections  are 

45  applied  ( Bondar  et  al,  1999).  Therefore,  energy  peaks  In  frequency-wavenumber 
44  spectra  observed  at  CMAR  can  be  confidently  backprojected  to  general  geo- 
e«  graphical  regions. 


51 


3  Computing  Frequency- Wavenumber  Spectra  at  CMAR 


k>  We  extracted  955  time  windows  of  seismic  noise  recorded  at  CMAR  from  1995 
u  through  2004.  The  160-s  long  windows  were  ch<:6en  to  end  several  tens  of  sec- 
»  onds  in  advance  of  P  waves  from  a  population  of  regional  -distance  earthquakes 
7®  that  were  analyzed  in  a  related  study.  The  time  windows  occur  in  all  months, 
n  years,  and  hours  of  the  day.  and  have  a  quasi-random  distribution.  All  the 
n  data  were  visually  inspected  and  channels  that  had  glitches  or  null  segments 
n  were  removed.  An  average  of  15.2  channels  per  noise  window  were  acceptable. 

74  We  used  two  complementary  methods  to  study  the  slowness  of  the  CMAR 
n  seismic  noise.  The  first  is  a  classic  frequency  domain  method  (Capon,  1969) 
n  that  averages  over  time  and  has  relatively  high  resolution  in  frequency.  The 
n  second  is  a  time  domain  beampacking  method  (e.g.,  Schweitzer  et  al.,  2002) 
it  that  averages  over  frequency  and  has  relatively  high  resolution  in  time.  For 
n  both  methods  we  assumed  that  the  seismic  noise  was  stationary  in  time  and 
»  space  over  the  sample  interval.  We  also  ignored  elevation  differences  among  the 
u  an  ay  elements  and  considered  only  2D  wavenumber  vectors.  This  is  appropri- 
k  ate  because  CMAR  is  relatively  flat,  having  a  maximum  elevation  difference  of 
m  only  0.066  km,  compared  to  a  horizontal  aperture  of  10. 1  km  (see  Bokelmann 
m  (1995)  for  a  discussion  of  topographic  effects  on  slowness  inference). 

«e  The  frequency  domain  technique  is  a  high-reeolution,  direct  segment  approach 
w  (Capon,  1969)  described  as  follows.  For  each  of  the  18  channels  of  CMAR  data 
w  we  extract  3200  data  points  and  subdivide  them  into  25  non-overlapping, 
»  adjacent  subwindows  of  128  samples  each.  The  sampling  interval  at.  CMAR  is 
»  0.05  s,  so  each  subwindow  is  6.4  s  long  and  the  total  window'  length  is  160  s.  For 
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»  each  subwindow  and  channel  the  mean  is  removed,  a  Bartlett  taper  is  applied, 
n  and  a  Fast  Fourier  Transform  (FFT)  is  calculated.  For  each  subwindow,  the 
k  18-by-18  spectral  matrix  is  calculated  by  multiplying  the  FFT  of  each  channel 
n  by  the  complex  conjugate  of  the  FFT  of  every  other  channel.  These  25  spectral 
«  matrix  estimates  are  then  linearly  averaged  to  produce  a  single  estimate  with 
«  reduced  \sriance.  Also,  by  averaging  more  subwindows  (25)  than  channels  (18) 
k  the  spectral  matrix  estimate  will  normally  become  nonsingular.  Finally,  the 
sc  matrix  is  normalized  by  dividing  each  element,  Sj{,  by  sJSj}Su.  This  corrects 
<e  for  unequal  gain  among  the  seismometers  at  CMAR. 

«  For  a  given  wavenumber  vector,  k  =  (kx,kv),  and  frequency,  fp,  the  power 
3oo  spectrum  is  defined  as, 

-  PUf  k)  =  (1) 


itc  where  Xj  is  the  position  vector  of  the  j,h  seismometer  relative  to  the  array 
3i»  reference  point,  and  J  is  the  number  of  channels  (18).  This  expression  can 
304  also  be  written  as  a  single  sum  if  the  set  of  all  spatial  difference  vectors, 
3k  called  the  coarray,  is  defined.  The  element  weights,  w},  define  a  2D  spatial 
sc*  windowing  function  that  is  analogous  to  the  Bartlett  taper  applied  in  the  time- 
307  domain  before  the  FFT.  Here  we  use  expressions  developed  for  a  maximum 
300  likelihood  filter  that  passes  energy  at  a  given  wavenumber  in  an  optimal  least- 
300  squares  sense  (Capon,  1969).  Discussions  of  various  strategies  for  chocsing 
330  these  weights  are  given  elsewhere  (Burg,  1964;  Haubrich.  1968;  Ingate  et  al., 
Ill  1985;  Douglas,  1998). 

332  With  the  time  windows  chcsen  in  this  work,  the  frequency  interval  at  which 
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us  P(fp,  k)  is  calculated  is  0.16  Hz.  For  a  given  frequency,  we  define  a  wavenum- 
n«  ber  grid  such  that  the  corresponding  slowness  grid  varies  from  -40  s/deg  to 
ns  40  s/deg  in  increments  of  0.5  s/deg,  and  evaluate  P(fp,  k)  and  at  each  node. 
ii»  The  grid  is  set  up  in  Cartesian  coordinates  and  aligned  such  that  the  x  compo- 
lw  nent  of  slowness  (sx  =  kx/  fp)  is  East/West,  and  the  y  component  of  slowness 
i»  (s„  =  ky/  fp)  is  North/South.  The  backazimuth  is  defined  clockwise  from  North 
i»  as  atan(Sj./Sy),  and  the  raj'  parameter  is  defined  in  units  of  s/deg  as  4-  sj. 
i»  We  also  define  the  apparent  or  phase  velocity  in  units  of  km /s  as  the  inverse 
22a  of  the  ray  parameter. 

us  The  time  domain  technique  is  based  on  algorithms  previously  developed  for 
us  analysis  of  coda  waves  and  small  amplitude  body  waves  (Koper  et  al,  2004). 
i3«  In  this  approach  we  first  extract  time  window's  of  length  64  s  from  each  of 
i3s  the  18  channels  at  C’MAR.  For  each  trace  wre  remove  the  mean,  apply  a  3 
i3»  pole  Butterworth  bandpass  filter,  and  resample  the  data  from  0.05  s  to  0.01  s. 
is?  The  comer  frequencies  of  the  filter  are  either  0.33-0.67  Hz.  0.67-1.33  Hz,  or 
i3*  1.33-2.67  Hz.  A  skwness  grid  with  the  same  dimensions  as  used  previously  is 
is*  constructed  and  a  phase  stack  weighted  beam  (Schimmel  <fc  Paulssen,  1997)  of 
i3»  order  three  is  formed  at  each  node.  The  power  is  defined  as  the  mean-square 
i3i  amplitude  of  the  beam  over  the  64  s  time  window.  This  approach  is  more 

133  computationally  intensive  than  the  frequency  domain  approach,  which  is  wrhy 
is3  we  use  a  shorter  time  window  (64  s  vs.  160  s).  In  all  cases,  the  window  for 

134  the  time  domain  method  is  chc6en  as  the  first  64  s  of  the  corresponding  160  s 
i3*  window. 

2se 

i3i  Example  spectra  axe  shown  in  Figure  2  for  a  typical  time  window'.  Power  is 
i3»  plotted  in  decibels  relative  to  the  sample  maximum,  10  x  logi0(P f i™,*).  and 
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i»  the  color  scale  is  chceen  such  that  only  energy  at  least  10%  of  the  maximum  is 
i«  highlighted.  Agreement  between  the  techniques  is  good,  especially  considering 
i«  that  the  time  and  frequency  sensitivities  are  different.  The  time-domain  tech- 
14:  nique  gives  slightly  sharper  images,  probably  because  it  amplifies  coherent 
i«  energy  via  phase-weighted  stacking.  In  both  methods,  at  higher  frequencies 
I**  the  noise  power  becomes  diffuse  and  unorganized  while  at  lower  frequencies 
i«e  there  are  specific  preferred  directions.  As  we  discuss  later,  the  preferred  noise 
i*«  directions  at  CMAR  are  generally  stable  from  one  time  window  to  the  next, 
i*j  though  the  absolute  amplitudes  exhibit  strong  variability. 


i*e  4  Ocean  Wave  Heights  From  TOPEX/POSEIDON  Data 

I**  To  help  interpret  noise  peaks  observed  in  fk-spectra  at  CMAR  we  used  a  global 
i»  model  of  ocean  wave  height  determined  from  satellite  data.  Since  ocean  waves 
m  play  a  dominant  role  in  generating  micrceeismic  noise,  it  follows  that  wave 
isi  heights  in  oceanic  source  regions  should  have  high  correlation  with  time  series 
i»  of  the  corresponding  microseismic  energy.  On  the  other  hand,  oceanic  regions 
is*  with  wave  heights  that  show  poor  or  negative  correlations  with  microseismic 
is*  energy  can  be  eliminated  as  potential  source  regions. 

is*  Significant  Wave  Heights,  defined  as  the  mean  height  of  the  highest  33%.  of 
is?  ocean  waves,  were  obtained  from  the  Ku-band  of  the  TOPEX/POSEIDON 
is*  satellite  altimetry  mission  (Callahan  et  al„  1994)  via  the  merged  geophysical 
]••  data  record  (MGDR),  distributed  by  the  NASA  Jet  Propulsion  Laboratory 
ie*  (JPL).  The  satellite  has  a  10-day  repeat,  orbit  with  a  maximum  orbiting  lati- 
i*i  tude  of  66°,  and  data  reported  every  second  are  averages  of  the  values  obt  ained 
in  from  10  Hz  measurements.  Using  data  for  which  at  least.  6  values  were  used 
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i a  in  the  average,  and  for  which  the  root  mean  square  error  of  the  6  or  more 
i«4  value6  was  below  1  m,  we  calculated  average  monthly  ocean  wave  heights  over 
iee  a  2°  by  2°  grid  for  the  same  ten  year  period  (1995-2004)  as  the  seismic  noise 
ie»  data.  Excluding  grids  that  contained  some  fraction  of  land  cover,  and  those 
it-  that  were  beyond  ±60:'  in  latitude  and  90  less  robust,  this  led  to  8,488  distinct 
iee  time  series  of  ocean  wave  height. 


ie»  5  Interpretation  of  FK-spectra  at  CMAR 

no  5.1  Average  Results 

m  Noise  power  as  a  function  of  slowness  is  shown  in  Figure  3  as  an  average 
m  over  the  955  time  windows.  Individual  spectra  were  normalized  before  con- 
in  tributing  to  the  stack,  therefore  the  relatively  few  windows  that  inadvertently 
n«  contain  large,  earthquake  generated  signals  do  not  bias  the  overall  result.  At 
in  frequencies  above  about  1.4  Hz  the  noise  is  incoherent  and  isotropic,  and  the 
m  two  methods  give  dissimilar  results.  However,  at  lower  frequencies  the  seis¬ 
in  mic  noise  is  strongly  anisotropic,  relative  to  both  backazimuth  and  apparent 
m  velocity,  implying  that  it  is  organized  as  propagating  plane  waves. 

m  At  lower  frequencies  the  two  methods  give  similar  results,  showing  three  dis- 
i»  tinct,  robust  peaks.  These  peaks  are  broadened  somewhat,  because  each  noise 
in  source  varies  within  a  general  geographic  region,  however  each  of  the  three 
ik  noise  sources  is  clearly  localized  as  shown  by  the  sharpness  of  the  correspond- 
iii  ing  histograms  of  peak  noise  power  (Figure  4).  In  that  figure  we  also  show  the 
im  effect  of  the  slowness  corrections  of  Bondar  et  al.  ( 1999).  We  implemented  the 
ik  default  correction  to  all  observations  and  further  applied  the  directionally  de- 
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I*  pendent  corrections  to  those  observations  falling  within  the  appropriate  bins 
>e?  in  slowness  space.  As  expected,  the  corrections  have  little  effect  on  the  overall 
i»  pattern  of  noise  peaks. 

)t*  It  is  important  to  emphasize  that  none  of  the  three  noise  peaks  shewn  in  Fig- 
)w  ures  3  and  4  are  sidelobes  of  one  another.  This  is  evident  by  examining  the 
m  CMAR  array  transfer  function,  shown  in  Figure  lb  at  1  Hz.  and  realizing  that 
i*2  at  lower  frequencies  the  sidelobe6  are  even  more  distant  in  slowness  space, 
iw  Furthermore,  both  of  the  spectral  estimation  techniques  we  used  reduce  the 
ih  amplitude  of  spurious  peaks,  with  the  frequency  domain  technique  in  partic¬ 
le  ular  essentially  deconvolving  the  array  transfer  function  from  the  observed 
i«t  spectra  (Capon,  1969). 

ik-  5.2  Rayleigh  Noise  from  the  East 

im  The  weakest  of  the  three  noise  peaks  shewn  in  Figure  3  has  apparent  velocities 
i»  centered  near  4.0  km/s  (indicative  of  higher  mode  Rayleigh  wraves)  and  arrives 
from  the  east  at  backazimuths  of  SO5  -  120°.  The  seismically  estimated  source 
2*i  area  for  this  Rayleigh  energy  is  shown  in  Figure  5  as  a  w'edge-shaped  region 
2*2  extending  eastward  from  CMAR  across  Vietnam  into  the  South  China  Sea. 
2*5  We  arbitrarily  bounded  the  source  region  at  a  distance  of  20"  because  of  the 
2«  high  attenuation  of  short-period  Rayleigh  w'aves. 

2*e  Surface  waves  are  commonly  observed  in  microseismic  noise  and  excitation  of 
2»  the  double  frequency  mode,  centered  at  periods  of  6-10  s,  can  be  explained 
2i*  by  counter-propagating  ocean  wraves  that,  interfere  in  a  nonlinear  manner 
2»  (Longuet-Higgins,  1950).  Our  observations  fall  on  the  short-period  side  of 


57 


this  peak,  but  can  nevertheless  be  explained,  for  instance,  by  the  interference 
no  of  ocean  waves  incident  on  and  reflected  by  the  Vietnamese  coastline.  Little 
m  seismic  noise  is  observed  from  the  northern  Vietnamese  coast  in  the  Gulf  of 
m  Tonkin,  which  is  presumably  shielded  from  the  Pacific  Ocean  by  the  island  of 
no  Hainan.  Likewise,  the  main  island  of  Malaysia  probably  acts  as  a  buffer  to  the 
no  south. 

ns  The  time  dependence  of  the  seismic  noise  power  is  shown  by  the  black  curve 
no  in  the  bottom  of  Figure  5.  For  each  time  window'  we  calculated  the  maximum 
n?  noise  power  in  the  relevant  region  of  slowness  space  and  then  normalized  each 
ne  estimate  by  the  median  of  the  population.  We  removed  nearly  20  samples  that 
n»  were  more  than  20  dBs  above  the  median,  to  eliminate  the  effect  of  earthquake 
n*  energy,  and  averaged  the  remaining  samples  into  monthly  bins.  The  signal  is 
m  strongly  seasonal  with  peaks  in  local  winter  and  troughs  in  local  summer,  and 
ns  has  a  dynamic  range  of  10-15  dB.  As  shown  in  Figure  5,  the  seasonality  of 
ns  the  seismic  noise  power  is  w'ell  matched  by  seasonal  variations  in  ocean  wrave 
n*  height  in  the  South  China  Sea.  The  sample  correlation  coefficient  betw'een  the 
ns  two  time  series  is  0.58,  which  is  above  the  99.9%  level  of  significance  using  a 
n*  Student-t  te6t. 

n?  5.3  Rayleigh  Noise  from  the  Soutkwest 

n«  The  largest  noise  peak  shown  in  Figures  3  arrives  from  the  southwest  at  back- 
ne  azimuths  of  ISO5  -  225°  and  also  has  phase  velocities  centered  near  4.0  km/s. 
n«  The  inferred  source  region  is  indicated  by  the  wedge  shape  in  Figure  6.  Ge- 
m  Graphically  it  corresponds  to  the  Burmese  coastline  from  the  border  with 
m  Thailand  in  the  south  to  the  peninsular  Pagoda  Point  in  the  north.  At  slightly 
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as  larger  baokazimuths  of  230°  -  280°  relatively  little  seismic  noise  is  generated 
2M  even  though  the  nearby  coastlines  of  western  Burma,  Bangladesh,  and  north- 
2x  eastern  India  lie  in  this  direction.  Likewise,  there  is  relatively  little  seismic 
a*  noise  generated  at  baokazimuths  of  135°  —  180°,  corresponding  to  coastlines 
a?  of  the  Malay  Peninsula  and  the  Gulf  of  Thailand.  These  two  seismically  quiet 
ae  regions  are  likely  shielded  from  northeasterly  travelling  storms  by  Sri  Lanka 
a»  and  Indonesia  respectively,  while  the  naiTow  region  in  between  is  subjected 
24c  directly  to  the  prevailing  currents  and  swell  of  the  larger  Indian  Ocean. 

s«  The  time  dependence  of  the  Indian  Ocean  seismic  noise  was  calculated  in  the 

342  same  manner  as  for  the  Pacific  Ocean  noise  and  is  shown  at  the  bottom  of 

343  Figure  6.  The  median  of  these  values  is  4.1  times  larger  than  the  median  of 
244  the  Pacific  Ocean  values.  Two  factors  are  likely  responsible  for  the  increased 
2«  power.  First,  there  is  larger  fetch  across  the  Indian  Ocean  than  the  South 
24t  China  Sea  so  the  excitation  of  the  Rayleigh  energy  on  the  southern  Burmese 

coast  is  probably  stronger  than  the  excitation  along  the  Vietnam  coast.  Sec- 
24t  ond,  the  seismic  propagation  path  from  coastline  to  CMAR  is  shorter  for 
24»  the  Indian  Ocean  noise  source,  and  so  lower  attenuation  of  the  short-period 
at  Rayleigh  energy  is  expected. 

ai  Like  the  Pacific  Ocean  noise,  the  Indian  Ocean  noise  is  strongly  seasonal  with 
as  a  dynamic  range  of  10-15  dB;  however,  the  seasonal  variation  is  opposite, 
as  with  peaks  in  local  summer  and  troughs  in  local  winter.  This  phenomenon 
a*  is  well  explained,  though,  by  the  seasonal  variation  of  wave  heights  in  the 
a*  Indian  Ocean,  which  peak  in  the  northern  summer  (southern  winter)  because 
a*  the  Indian  Ocean  as  a  wrhole  has  southern  hemisphere  dynamics.  Thus  the 
a?  seasonal  anticor relation  of  wrave  heights  in  the  Pacific  and  Indian  Oceans  is 
at  mimicked  by  the  seasonal  anticorrelation  of  seismic  noise  from  these  regions. 
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5.4  Tele  seismic  Body  Wave  Noise 

m®  The  third  noise  peak  shown  in  Figure  3  occurs  near  the  oenter  of  the  diagram, 
aw  indicating  energy  that  is  almost  vertically  incident  on  the  array.  Ray  param- 
xi  eters  for  this  peak  are  clustered  between  0  and  5  s/deg,  with  a  mean  near 
»3  2.5  s/deg  for  the  uneorrected  data  (Figure  4c)  and  a  mean  near  3.5  s/deg  for 
3®4  the  corrected  data  ( Figure  4d).  Owing  to  the  near- vertical  incidence  the  back- 
3«  azimuths  are  poorly  resolved,  though  there  is  a  trend  towards  the  southeast. 

The  ray  parameter  for  P  waves  that  turn  just  above  the  core-mantle  bound- 
2*7  ary  (equivalently  P<h> / )  is  about  4.4  s/deg,  so  the  fact  that  we  observe  smaller 
2M  values  implies  that  the  noise  is  composed  of  P  waves  that  have  interacted  with 
:»  Earth's  core. 

27®  We  first  consider  the  possibility  that  the  waves  have  been  transmitted  through 
27i  the  core  as  PKP  wave6,  and  assume  that  the  noise  source  occurs  at  a  distance 
371  of  140°  -  145°  from  CMAR,  as  shown  in  Figure  7.  This  range  brackets  the 
371  PKP  -  B  caustic,  a  distance  at  which  energy  is  geometrically  focused  to  a 
274  maximum  by  the  velocity  structure  in  Earth’s  core.  We  note  that  previous 
27!  'observations  of  P  waves  in  microseismic  noise  have  generally  occurred  at  raj' 
37®  parameters  consistent  with  a  distance  of  20°  (Haubrich  &  McCamy,  1969; 
377  Gerstoft  et  al„  2006),  a  range  in  which  body  wrave  energy  is  also  geometrically 
it®  focused,  in  this  case  by  high  velocity  gradients  in  the  mantle  transition  zone. 

it®  The  time  dependence  of  the  P-wrave  noise  is  shewn  by  the  black  line  in  the 
2*®  bottom  of  Figure  7.  It  wras  calculated,  in  the  manner  described  earlier,  from  the 
2®i  region  of  slowness  space  with  velocity  greater  than  20  km/s  (ray  parameter 
3*1  less  than  5.5  s/deg),  and  it  has  a  median  that  is  3.1  times  larger  than  the 
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»j  median  of  the  Pacific  Ocean  Rayleigh  noise.  There  is  a  clear  seasonality  to 
■m  the  seismic  noise,  with  peaks  in  the  northern  winter  and  troughs  in  northern 
:«  summer,  though  the  dynamic  range  is  only  5-10  dB.  The  correlation  of  the 
:»  seismic  noise  with  ocean  wave  heights  acrcss  the  hypothetical  PKP  source 
»-■  region  is  also  shown  in  Figure  7.  While  it  is  tempting  to  suggest  that  the 
:*«  strong  waves  near  the  southern  tip  of  South  America  create  PKP  noise  at 
■m  CMAR,  no  significant  correlation  was  found  with  the  seismic  observations. 
:«c  The  only  portion  of  the  seismic  source  annulus  that  gives  high  correlation 
»i  is  the  patch  along  the  northern  coast  of  Brazil  near  the  Amazon  delta.  The 
7K  sample  correlation  coefficient  for  this  area  is  0.44,  which  is  significant  above 
the  99.9%  level. 

*4  A  potential  problem  with  this  interpretation  though  is  that  it  corresponds  to 
■m  northwesterly  backazimuths  at  C’MAR  The  data  prefer  southeasterly  backaz- 
:*  imuths  ( Figure  4),  corresponding  to  the  southwestern  quadrant  of  the  seismic 
jw  source  annulus  shown  in  Figure  7.  Therefore,  for  this  interpretation  to  be  eor- 
:««  rect  it  must  be  the  case  that  our  slowness  estimates  are  biased  (on  the  order 
:*»  of  3-4  s/deg)  by  imperfections  in  the  slowness  corrections  of  Bondar  et  al. 
see  (1999),  or  because  of  some  deficiency  in  our  observational  techniques.  For  ex- 
*»  ample,  our  slowness  inference  techniques  assume  a  single  plane  wave  model 
a «  and  if  multiple  plane  wrave6  arrive  in  the  same  time  window  it  is  possible  that 
in  the  peak  of  the  dominant  arrival  is  shifted  somewhat  in  slowness  space  (e.g., 
304  Shunway  et  al.,  2008). 

3«  We  next  consider  the  possibility  that  the  high-velocity  noise  peak  observed  at 
so*  CMAR  is  created  by  energy  that  has  reflected  o-ff  the  core-mantle  boundary  as 
3w  PcP  waves.  This  seems  less  plausible  than  the  PKP  hypothesis  because  PcP 
so*  waves  are  relatively  small,  especially  at  small  source-receiver  distances,  and 
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x*  the  corresponding  direct  P  waves  would  be  expected  to  be  visible  in  the  noise 
no  spectra  with  raj’  parameters  of  6-9  s/deg.  Nevertheless,  it's  plausible  that,  at 
su  least  some  of  the  high-velocity  noise  observed  at  CMAR  is  created  in  this 
in  manner. 

ns  In  Figure  8  we  show'  the  correlation  of  the  high-velocity  noise  time  series  with 
314  ocean  wave  heights  for  distances  less  than  90°  from  C’MAR,  which  is  equivalent 
ns  to  the  hemisphere  opposite  that  shown  in  Figure  7.  The  seemingly  strong  eorre- 
3it  lation  in  the  Barents  Sea  is  an  artifact  related  to  the  extrapolation  of  our  ocean 
3i?  height  model  to  latitudes  above  what  is  covered  by  TOPEX/POSEIDON  satel- 
3i>  lite  tracks;  however,  the  strong  correlations  observed  across  much  of  the  north 
114  Pacific  are  legitimate.  In  the  bottom  panel  we  compare  the  seismic  noise  time 
s»  series  with  wave  heights  from  a  domain  just  north  of  New  Guinea,  at  a  dis- 
3n  tance  of  aro>und  50°  from  CMAR.  A  PcP  wave  from  this  area  would  arrive 
in  at  CMAR  with  a  ray  parameter  of  3.7  s/deg  from  a  backazimuth  of  around 
sis  100°.  The  sample  correlation  coefficient  in  this  case  is  0.37,  which  is  significant 
324  above  the  99.9%  level.  Equivalently  strong  positive  correlations  are  found  for 
ns  regions  in  the  North  Pacific  near  the  Aleutians,  a  region  known  for  large  wave 
Ht  heights. 

su  A  final  possibility  is  that  the  high-velocity  noise  recorded  at  CMAR  is  created 
32*  by  direct  teleseismic  P  waves  that  turn  deep  in  the  lower  mantle.  High  corre- 
i2»  lation  with  ocean  w'aves  in  the  Pacific  is  seen  for  many  regions  at  distances  of 
s»  60°  -  90:i  from  CMAR  (Figure  8).  These  P-wave6  would  arrive  at  CMAR  from 
in  the  east  with  ray  parameters  of  4. 5-7.0  s/deg,  contradicting  our  observations 
i32  by  about  3-4  s/deg  in  absolute  slowness.  But  this  is  the  same  margin  of  error 
su  as  for  the  PKP  interpretation,  and  so  for  the  reasons  mentioned  earlier  the 
314  teleseismic  P  energy  could  plausibly  contribute  to  the  observed  noise. 
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6  Conclusions 


J36 


j»  We  calculated  frequency-wavenumber  spectra  for  955  samples  of  seismic  noise 
m?  recorded  at  the  Chiang  Mai  seismic  array  (CMAR)  from  1995  through  2004. 
s»  At  frequencies  above  about  1.4  Hz  the  noise  is  isotropic  and  diffuse.  Therefore 
:»  in  array  analysis  of  body  waves  from  small  regional  events  recorded  at  CMAR, 
j**  the  assumption  of  an  isotropic  noise  model  is  justified  and  gains  close  to  %/l8 
mi  are  expected  when  processing  high-passed  waveforms  with  coherent  signals. 

mi  At  frequencies  lower  than  1.4  Hz  the  noise  at  CMAR  is  strongly  partitioned  by 
ms  apparent  velocity  into  two  categories:  teleeeismic  P  wave  energy'  with  apparent 
344  velocities  higher  than  25  km/s  (ray  parameters  of  0.0-5.0  s/deg)  and  higher 
s4s  mode  Rayleigh  energy  with  apparent  velocities  near  4.0  kra/s.  The  ring  of 
34*  slowness  space  in  between,  corresponding  to  P  waves  tinning  in  the  crust  and 
34?  upper  mantle,  is  relatively  quiet.  The  Rayleigh  noise  is  further  partitioned 
34®  by  direction,  with  the  strongest  signal  arriving  from  the  Bay  of  Bengal  at 
34.  backazimuths  of  1805  —  255°.  A  secondary  peak  in  the  Rayleigh  noise  occurs 
s»  in  the  direction  of  the  South  China  Sea  at  backazimuths  of  80°  -  120°.  Little 
sei  noise  arrives  from  the  north,  the  direction  of  the  main  Asian  landmass  and 
3d  the  nuclear  test  sites  of  India,  Pakistan.  China,  and  North  Korea. 

363  The  Rayleigh  noise  is  strongly  seasonal  with  annual  variations  of  10-15  dB 

364  in  power .  The  easterly  noise  has  peaks  in  local  winter  and  troughs  in  lo¬ 
se?  cal  summer,  while  the  noise  from  the  southwest  has  the  opposite  pattern. 
36«  This  behavior  is  well-matched  by  the  seasonal  anti-correlation  in  significant 
se?  wave  heights  in  the  South  China  Sea  and  the  Bay  of  Bengal,  as  determined 
36»  from  TOPEX/POSEIDON  satellite  tracks.  For  a  magnitude  scale  such  as 
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*«  MS(Vmax)  (Russell,  2006;  Bonner  et  al.,  2006)  that  depends  directly  on  the 
j»  logarithm  of  seismic  amplitude,  our  observations  imply  that  the  detection 
jw  threshold  at  CMAR  varies  b>'  0.5-0.75  magnitude  units  according  to  direction 
aes  and  time  of  year. 

3»j  While  propagating  Rayleigh  waves  are  often  observed  in  seismic  noise,  it  is 
364  less  common  to  observe  teleseismic  body  waves.  Nearly  all  reports  of  body 
see  wave  noise  document  ray  parameters  of  8-12  s/deg,  which  correspond  to  P 
3 «  waves  that  turn  in  the  upper  mantle  (e.g.,  Toksoz  &  Laccss,  1968;  Haubrich 
3*.-  &  McCamy,  1969;  Cole  et  al.,  1989;  Gerstoft  et  al,  2006).  Therefore,  our 
m  observations  of  a  consistent  noise  peak  with  ray  parameters  of  4.5  s/deg  and 
3»  smaller,  equivalent  to  apparent  velocities  of  25  km/s  and  higher,  may  be  unique 
m  in  the  geophysical  literature. 

in  Like  the  Rayleigh  noise,  the  P  noise  observed  at  CMAR  is  seasonal.  It  has 
372  an  annual  power  variation  of  5-10  dB,  with  peaks  in  local  winter  and  troughs 
vs  in  local  summer.  The  seasonality  implies  that  the  noise  is  unrelated  to  small 
374  unidentified  earthquakes  or  other  tectonic  processes,  and  instead  is  created 
37i  at  least  indirectly  by  ocean  waves,  similar  to  the  Rayleigh  noise.  Under  this 

376  assumption  there  are  several  geographical  regions  that  could  act  as  sources: 

377  the  western  Atlantic  Ocean  near  the  coast  of  northern  Brazil  may  contribute 
376  PKP  energy,  the  Pacific  Ocean  just  north  of  New  Guinea  may  contribute  PcP 
376  energy,  and  central  portions  of  the  North  Pacific  may  contribute  P  waves  that 
se*  turn  in  the  lower  mantle.  However,  none  of  these  sources  provides  an  ideal 
36i  match  to  our  slowness  observations  and  so  none  are  preferred  over  the  others. 

36i  In  order  to  resolve  the  uncertainty  of  the  source  location  of  the  P  noise  at 
36i  CMAR,  a  finer  comparison  of  seismic  data  and  ooean  wave  data  is  required. 
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Ideally,  we  would  compare  dally  or  even  hourly  samples  of  seismic  noise  to 
»e  the  geographical  tracks  of  named  storms  and/or  to  equally  uniformly  sampled 
;«  ocean  wave  height  data  from  buoys.  Observations  of  P  noise  from  other  arrays 
»  in  the  general  area  would  also  help,  especially  if  these  arrays  are  slightly  wider 
j«  aperture  than  CMAR  with  correspondingly  smaller  slowness  biases  related  to 
m»  near-receiver  heterogeneity.  In  any  case,  and  irrespective  of  the  precise  source 
s»  mechanism  for  the  high  velocity  noise,  our  observations  point  towards  a  new 
method  of  imaging  Earth's  deep  interior.  Just  as  Rayleigh  microeeismic  noise 
sk  has  been  used  to  image  Earth's  crust  (Shapiro  et  al,  2005)  it  may  be  possible 
3«  to  use  micrcseismic  P-noise  to  image  Earth's  lowermost  mantle  and  core. 
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:«  Data  and  Resources.  The  TOPEX/ POSEIDON  Significant  Wave  Height 
x  data  were  obtained  from  the  Physical  Oceanography  Distributed  Active  Archive 
*»  Center  (PO.DAAC)  at  the  NASA  Jet  Propulsion  Laboratory,  Pasadena,  CA. 
jw  They  ate  publically  available.  Seismograms  used  in  this  study  were  obtained 
we  from  U.S.  National  Data  Center  ( NDC)  maintained  by  the  Air  Force  Technical 
Me  Applications  Center  (AFTAC)  in  Melbourne,  Florida.  They  are  also  available 
4 tC  to  the  public. 
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Figure  Captions 

4«4 

«  Figure  1.  (left)  Location  and  geometry  of  the  C’hiang  Maj  seismic  array 
4w  (CMAR).  The  dark  triangle  is  the  reference  point  for  the  array,  (right)  The 
w  array  response  function  for  a  vertically  incident  1  Hz  plane  wave.  Note  that 
«e  the  sidelobes  are  distant  and  weak.  At  lower  frequencies,  the  sidelobe6  are 
«*  even  further  separated  from  the  main-lobe  in  slowness  space. 

49C 

«i  Figure  2.  Frequency- wavenumb er  spectra  for  a  time  window  starting  of  seis- 
4«  mic  noise  starting  at  approximately  08:30:00  GMT  on  March  26,  1906.  The 
panels  on  the  left  shew  results  from  the  frequency  domain  technique  for  a  160  s 
«*  time  window,  and  the  panels  on  the  right  shew  results  from  the  time  domain 
«  technique  for  a  64  s  subwindow.  Frequency  increases  frem  top  to  bottom  in 
w«  both  cases.  The  inner  ring  indicates  a  ray  parameter  of  4.4  s/deg  (phase  veloc- 
«•  ity  of  about  25  km/s),  the  middle  ring  indicates  a  raj'  parameter  of  27.8  s/deg 
w*  (phase  velocity  of  4  km/s),  and  the  outermost  ring  indicates  a  raj’  parameter 
4»  of  40.0  s/deg  (phase  velocity  of  about  2.8  km/s).  Therefore,  energy  arriving 
!»  within  the  inner  ring  possesses  a  ray  parameter  smaller  than  and  so 
sci  has  likely  interacted  with  Earth’s  core.  Energy  arriving  between  the  two  outer 
sos  rings  is  consistent  with  higher-mode  Rajdeigh  waves,  and  is  typical  of  short- 
scs  period  noise  recorded  at  ar raj's. 

!«4 

s«  Figure  3.  Stacked  frequency  wavenumber  spectra  for  the  955  time  windows, 
so*  Subpanels  are  the  same  as  described  in  Figure  2. 
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k«  Figure  4.  Histograms  of  peak  noise  power  for  the  955  time  windows.  The 
•»  global  maxima  were  selected  from  the  time-domain  algorithm  applied  in  the 
mo  0.33-0.67  Hz  band.  We  show'  the  full  2D  histogram  in  0.5  s/deg  by  0.5  s/deg 
on  bins  with  circles  drawn  at  apparent  velocities  of  25  km/s,  4  km/s,  and  2.8  km/s 
as  for  the  (a)  raw  and  (b)  corrected  slownesses.  The  corresponding  raj’  parame- 
5is  ter  histogram  in  bins  of  1  s/deg  are  shown  in  panels  (c)  raw  and  (d)  corrected. 

U4 

m  Figure  5.  Correlation  of  the  Pacific  Ocean  Rayleigh  wave  micrcseismic  energy 
ut  with  ocean  wave  heights.  In  the  top  panel  the  seismically  constrained  source 
a?  area  is  shown  by  the  wedge  shape,  and  the  color  indicates  the  correlation  co¬ 
ne  efficient  between  the  seismic  noi9e  and  the  significant  ocean  wave  heights.  The 
He  time  dependence  of  the  seismic  noise  is  shown  by  the  black  curve  in  the  bot¬ 
es®  tom  panel,  and  the  red  curve  represents  significant  wave  height  for  a  domain 
5si  of  11°N-15°N  and  111CE-115°E. 

822 

ess  Figure  6.  Correlation  of  the  Indian  Ocean  Rayleigh  wave  micrcseismic  energy 
53<  with  ocean  wave  heights.  In  the  top  panel  the  seismically  constrained  source 
535  area  is  shown  by  the  wedge  shape,  and  the  color  indicates  the  correlation  co- 
5s»  efficient  between  the  seismic  noise  and  the  significant  ocean  wave  heights.  The 
S3?  time  dependence  of  the  seismic  nc-ise  is  shown  by  the  black  curve  in  the  bot- 
53*  tom  panel,  and  the  red  curve  represents  significant  wave  height  for  a  domain 
53.  of  8’N-12’N  and  91°E-95°E. 

530 

55i  Figure  7.  Correlation  of  the  high-velocity  micrcseismic  energy  with  ocean 
si!  wave  heights,  assuming  a  PKP  source.  In  the  top  panel  the  seismically  con- 
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sss  strained  source  area  is  shown  by  the  annulus,  and  the  color  indicates  the  eor- 
ss4  relation  coefficient  between  the  seismic  noise  and  the  significant  ocean  wrave 
ms  heights.  The  time  dependenoe  of  the  seismic  noise  is  shown  by  the  black  curve 
Me  in  the  bottom  panel,  and  the  red  curve  represents  significant  wave  height  for 
mj  a  domain  of  near  the  north  coast  of  Brazil  defined  as  2:>S-3':’N  and  50°W-45°W. 

536 


t»  Figure  8.  Correlation  of  the  high-velocity  micrcseismic  energy'  with  ocean 
sec  wave  heights,  assuming  a  P/PcP  source.  In  the  top  panel  the  circles  are  drawm 

543  at  distances  of  30°  and  60,:  from  CMAR,  and  the  color  indicates  the  correlation 
542  coefficient  between  the  seismic  noise  and  the  significant  ocean  wave  heights. 
545  The  time  dependence  of  the  seismic  noise  is  showm  by  the  black  curve  in 

544  the  bottom  panel,  and  the  red  curve  represents  significant  wave  height  for 
54s  a  domain  just  to  the  north  coast  of  New  Guinea  defined  as  0CN-10CN  and 
54*  135°E-150CE. 
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